Can we approximate harmonic functions with harmonic functions with non-vanishing differential? Let $\mathbb{D}^2$ be the closed $2$-dimensional unit disk, and let $g:\mathbb{D}^2 \to \mathbb{R}$ be a non-constant harmonic function (in particular smooth up to the boundary).

Does there exist a sequence of smooth harmonic functions $g_n$ on $\mathbb{D}^2 $, such that $g_n \to g$ in $W^{1,2}$ and $dg_n \neq 0$ everywhere on $ \text{int}(\mathbb{D}^2)$?

Since we can add additive constants to the $g_n$, we can arrange $\int_{\mathbb D^2} g_n=\int_{\mathbb D^2} g$, so the $W^{1,2}$ convergence of the $g_n$ is essentially equivalent to $dg_n \to dg$ in $L^2$. (via Poincare inequality). 
Thinking on $dg$ as a vector field, I think that we can always approximate it with a non-zero vector field in $L^2$. However, the only procedure I know for doing that does not produce approximating vector fields which are gradients of harmonic functions (or gradients of anything, really).
 A: No, even if we only ask $g_n\to g$ in $L^2$ (or even just in the sense of distributions, but that gets more technical).
The important property is that $L^2$ convergence of harmonic functions on an open set implies uniform convergence of each derivative on each compact subset. This is a kind of topological version of Weyl's lemma: not only are harmonic functions smooth, but "convergent sequences of harmonic functions" are "convergent sequences of smooth functions" in a certain sense.
Set $g(x,y)=x^2-y^2.$
The mean value property of harmonic functions implies that $\int f(z'-z)\psi(z)dx=f(z')\int \psi$ whenever $f$ is harmonic on a ball $B(z',r),$ and $\psi$ is radially symmetric with support in the ball $B(0,r).$
Take a radially symmetric smooth function $\psi$ supported in the ball $B(0,1/2)$ and with $\int \psi = 1.$ For any $z\in B(0,1/2)$ we get
$$dg(z) = \int dg(z'-z)\psi(z) = -\int g(z'-z)d\psi(z)dz$$
and similarly for $g_n.$ So
$$|(dg-dg_n)(z)|\leq |\int (g(z'-z)-g_n(z'-z))d\psi(z)dz|\leq \|g-g_n\|_2\|d\psi\|_2$$
which tends to zero uniformly in $z$ as $n\to\infty.$
$dg$ has winding number $-1$ around $0$ as $z$ goes anticlockwise around along the circle $|z|=1/2.$ So for large enough $n,$ by continuity $dg_n$ gets the same winding number and therefore must have a zero in the disc $|z|\leq 1/2.$
(I think a similar idea should answer https://mathoverflow.net/questions/330947/can-we-approximate-harmonic-maps-which-are-a-e-orientation-preserving-with-maps: take $f(x,y)=(x^2-y^2,2xy)$ and apply the above argument to $f_1(x,y)=x^2-y^2.$)
