This works if and only if $\partial_j g_k=\partial_k g_j$ in distributional sense for all $j,k=1,\ldots , n$.

Obviously, since $g_k=\partial_k f$ in $\mathcal D'$, this condition is necessary. Conversely, if these integrability conditions hold, we can approximate in $L^p$ by smooth functions $h_n\to g$ that still satisfy them (for example, take convolutions with suitable smooth functions). By Fubini, we then also have $h_n(\cdot, y)\to g(\cdot,y)$ in $L^p$ for almost every $y$. Let's assume, for convenience, that $y=0$ works here, and, also for convenience, I'll assume $n=2$ from now on.

Let $f_n(x)=\int_{\gamma(x,y)} h_n(t)\, dt$ and $f(x)=\int_{\gamma(x,y)} g(t)\, dt$, with $\gamma$ denoting a path that goes from $(0,0)$ to $(x,y)$ along straight line segments parallel to the coordinate axes. Then \begin{align*} \|f_n-f\|_p^p & \le \int\!\int dx\, dy \left( \int_{\gamma(x,y)}dt\, |h_n(t)-g(t)|\right)^p \\ &\le \int\int dx\,dy\, (|x|+|y|)^{p-1} \int_{\gamma(x,y)} dt\, |h_n(t)-g((t)|^p\\ & = \int\int dx\,dy\, (|x|+|y|)^{p-1} \\ & \quad \times \left( \int dt\,|h_n(t,0)-g(t,0)|^p +\int dt\, |h_n(x,t)-g(x,t)|^p \right) . \end{align*} Here I have used Jensen's inequality and then we can compare the last double integral with how we would compute the norm in spherical coordinates.

It follows that $f_n\to f$ in $L^p$ also. (Strictly speaking, we should first run this calculation without $h_n$ to confirm that $f\in L^p$.) Moreover, $\nabla f_n= h_n\to g$ in $L^p$, so $f\in W^{1,p}$ and $\nabla f=g$, as desired.

This works if and only if $\partial_j g_k=\partial_k g_j$ in distributional sense for all $j,k=1,\ldots , n$.

Obviously, since $g_k=\partial_k f$ in $\mathcal D'$, this condition is necessary. Conversely, if these integrability conditions hold, we can approximate in $L^p$ by smooth functions $h_n\to g$ that still satisfy them (for example, take convolutions with suitable smooth functions). By Fubini, we then also have $h_n(\cdot, y)\to g(\cdot,y)$ in $L^p$ for almost every $y$. Let's assume, for convenience, that $y=0$ works here, and, also for convenience, I'll assume $n=2$ from now on.

Let $f_n(x)=\int_{\gamma(x,y)} h_n(t)\, dt$ and $f(x)=\int_{\gamma(x,y)} g(t)\, dt$, with $\gamma$ denoting a path that goes from $(0,0)$ to $(x,y)$ along straight line segments parallel to the coordinate axes. Then \begin{align*} \|f_n-f\|_p^p & \le \int\!\int dx\, dy \left( \int_{\gamma(x,y)}dt\, |h_n(t)-g(t)|\right)^p \\ &\le \int\int dx\,dy\, (|x|+|y|)^{p-1} \int_{\gamma(x,y)} dt\, |h_n(t)-g((t)|^p\\ & = \int\int dx\,dy\, (|x|+|y|)^{p-1} \\ & \quad \times \left( \int dt\,|h_n(t,0)-g(t,0)|^p +\int dt\, |h_n(x,t)-g(x,t)|^p \right) . \end{align*} Here I have used Jensen's inequality to pass to the second line.

It follows that $f_n\to f$ in $L^p$ also. (Strictly speaking, we should first run this calculation without $h_n$ to confirm that $f\in L^p$.) Moreover, $\nabla f_n= h_n\to g$ in $L^p$, so $f\in W^{1,p}$ and $\nabla f=g$, as desired.