Here is a short, elementary, and self-contained proof of the result you wanted to prove. It is similar to the one given in Maz'ya's book, but simpler.
For a related post see: https://mathoverflow.net/a/297392/121665
If $f\in L^1_{\rm loc}$ or even if $f$ is a distribution and $\nabla
f\in L^p$, then $f\in L^p_{\rm loc}$.
This follows from the fact proven below. For $1\leq p\leq \infty$ define
$$
L^{1,p}(\Omega)=\{ f\in \mathcal{D}'(\Omega):\, \nabla f\in L^p(\Omega)\},
\quad
W^{1,p}(\Omega)=\{ f\in L^p(\Omega):\, \nabla f\in L^p(\Omega)\}.
$$
Similarly we define spaces $L^{1,p}_{\rm loc}$ and $W^{1,p}_{\rm loc}$.
Theorem. $L_{\rm loc}^{1,p}(\Omega)\subset W^{1,p}_{\rm loc}(\Omega)$ for $1\leq p\leq\infty$.
Remark. By induction the result generalizes to spaces with higher order derivatives.
In the proof we will need the following well known result.
Lemma. If $u\in W^{1,1}(\mathbb{R}^{n})$ then $$ u(x)= \frac{1}{n\omega_{n}}\int_{\mathbb{R}^{n}}\frac{(x-y)\cdot \nabla u(y)}{|x-y|^{n}}\, dy \quad \text{a.e.,} $$ where $\omega_n$ denotes the volume of the unit ball.
Proof. By a density argument it suffices to prove it for $u\in C_0^\infty(\mathbb{R}^n)$. Let
$s\in S^{n-1}$ (unit sphere). We have
$$
u(x)=-\int_{0}^{\infty}\frac{d}{dr}u(x+rs)\, dr=
-\int_{0}^{\infty}Du(x+rs)\cdot s\, dr.
$$
Taking the average with respect to $s\in S^{n-1}$ we get
(recall that the volume of $S^{n-1}$ equals $n\omega_n$)
$$
u(x)=-\frac{1}{n\omega_{n}}\int_{S^{n-1}}\int_{0}^{\infty}Du(x+rs)\cdot
s\, dr\, ds
$$
and the lemma follows after substituting
$x+rs=y$, so $dy=r^{n-1}drds$, $s=(y-x)/|y-x|$,
and $drds=|x-y|^{(1-n)}dy$. $\Box$
The lemma can be seen as the
integral representation of the Dirac $\delta$ distribution.
Corollary.
$\sum_{i=1}^{n}\partial K_{i}/\partial x_{i}=\delta$, where $K_{i}=n^{-1}\omega_{n}^{-1}x_{i}|x|^{-n}$.
$\Box$
Proof of the theorem.
Suppose that $u\in L^{1,p}(\Omega)$. Let
$V\Subset V_{\varepsilon}\Subset\Omega$, where
$V_{\varepsilon}=\{ x\, |\, {\rm dist}\,(x,V)<\varepsilon\}$ and let
$\varphi\in C_0^\infty(\Omega)$, $\varphi|_{V_{\varepsilon}}\equiv 1$. It suffices to
prove that $w=\varphi u\in L^{p}(V)$. In other words it suffices to prove
that the distribution $\varphi u$ when restricted to $V$
can be represented by a certain $L^{p}(V)$ function.
Let $\eta\in C_0^\infty(B^{n}(0,\varepsilon))$,
$\eta|_{B^{n}(0,\varepsilon/2)}\equiv 1$. According to the corollary we have
$$
\sum_{i=1}^{n}\frac{\partial(\eta K_{i})}{\partial x_{i}}=
\sum_{i=1}^{n}\eta\frac{\partial K_{i}}{\partial x_{i}}+\xi=\delta+\xi,
$$
where
$\xi\in C_0^\infty(B^{n}(0,\varepsilon))$.
The fact that $\xi$ is smooth follows from the observation
$K_{i}\partial\eta/\partial x_{i}\in C_0^\infty(B^{n}(0,\varepsilon))$.
Now
$$
w+\xi*w=(\delta+\xi)*w=\sum_{i=1}^{n}\frac{\partial(\eta K_{i})}{\partial x_{i}}*w
=\sum_{i=1}^{n}(\eta K_{i})*\frac{\partial w}{\partial x_{i}}.
$$
By properties of the convolution, the distribution $\xi*w$ is
smooth in $\mathbb{R}^{n}$ and hence it belongs to $L^{p}(V)$ (after
being restricted to $V$). Therefore it remains to show that
$$
(\eta K_{i})*\frac{\partial w}{\partial x_{i}}\in L^{p}(V)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (*)
$$
Since ${\rm supp}\,\eta K_{i}\subset B^{n}(0,\varepsilon)$, the
restriction of the distribution
$(\eta K_{i})*(\partial w/\partial x_{i})$ to $V$ does not depend on the
behavior of $\partial w/\partial x_{i}$ outside $V_{\varepsilon}$, but
$\partial w/\partial x_{i}=\partial u/\partial x_{i}$ in
$\mathcal{D}'(V_{\varepsilon})$ (as $\varphi=1$ in $V_{\varepsilon}$) so
$$
(\eta K_{i})*\frac{\partial w}{\partial x_{i}}=
(\eta K_{i})*\frac{\partial u}{\partial x_{i}}
\quad
\text{in $\mathcal{D}'(V)$.}
$$
Now (*) follows since
$$
\Vert(\eta K_i)*\frac{\partial u}{\partial x_i}\Vert_p\leq\Vert\eta K_i\Vert_1
\Vert\frac{\partial u}{\partial x_i}\Vert_p\leq C \Vert\frac{\partial u}{\partial x_i}\Vert_p.
$$
Indeed, $\eta$ has compact support and $|K_i|\leq C|x|^{1-n}$ so $\eta K_i\in L^1$.
$\Box$
