I believe that Daniele Tampieri's idea of using completeness Fichera's Theorem is the right one. However I think there is a simpler proof.
In view of the Hahn-Banach Theorem,  proving the density of the space
\begin{equation}
\left\{(f,f|_{\partial D})\ |\ f\in S \right\}, \quad \text{where } 
S=\left\{ f\in C^{\infty}(\overline{D})\ |\ \Delta f|_{\partial D} =0 \right\},
\end{equation}
in $L^{2}(D)\times L^{2}(\partial D)$ 
is equivalent to show that, if $(F,G)\in L^{2}(D)\times 
L^{2}(\partial D)$ is such that
\begin{equation}
\int_{D}F\, f\, dx + \int_{\partial D}G\, f\, d\sigma = 0
\end{equation}
for any $f\in S$, then $F=0$ a.e. in $D$
and $G=0$ a.e. on $\partial D$.
So suppose that these orthogonality conditions hold. In particular we have
$$
\int_{D}F\, f\, dx   = 0
$$
for any $f\in C^{\infty}_{c}(D)$. It is well known that this implies
$F=0$ a.e. in $D$. 
Therefore
$$
\int_{\partial D}G\, f\, d\sigma = 0
$$
for any $f\in S$. Harmonic polynomials obviously belong to $S$ and then
$$
\int_{\partial D}G\, \omega\, d\sigma = 0
$$
for any harmonic polynomial $\omega$. Fichera's theorem implies
 $G=0$ a.e. on $\partial D$.