Gauss–Ostrogradsky formula for Distributions
Solution 1:
After a long search, I found the answer in Constantin's book. I will post a more general answer here. Let $\Omega\subset\mathbb{R}^N$ be a bounded domain with Lipschitz boundary. Let $\gamma_{0,p}:W^{1,p}(\Omega)\to W^{1-1/p,p}(\partial\Omega)$ be the trace operator and $\ell_p:W^{1-1/p,p}(\partial\Omega)\to W^{1,p}(\Omega)$ the lift operator.
The lift operator is continuous and determined by, $u=\ell_p v$ $$ \left\{ \begin{array}{ccc} \Delta_pu=0 &\mbox{ in $\Omega$} \\ \gamma_{0,p}u=v &\mbox{ in $\partial\Omega$} \end{array} \right. $$
Let $E_p=\{u\in L^p(\Omega)^N:\ \operatorname{div}u\in L^p(\Omega)\}$.
Note that $E_p$ is a Banach space with the norm $$\|u\|_E=\|u\|_p+\|\operatorname{div}u\|_p$$
Define $\gamma:E_p\to W^{-1+1/p,p}(\partial\Omega)^N$ by $$\langle \gamma (u),v\rangle_{W^{-1+1/p,p},W^{1-1/p',p'}}=\int_\Omega \ell_p(v)\operatorname{div}u+\int_\Omega \nabla(\ell_p (v))\cdot u $$
It follow that (by using the continuity of $\ell_p$) $$|\langle\gamma(u),v\rangle|\leq c\|u\|_E\|v\|_{W^{1-1/p',p'}} $$
which implies that $\gamma$ is a bounded linear operator. Now, we note that if $u\in C^\infty(\overline{\Omega})^N$ and $v\in C^\infty(\overline{\Omega})$, then, by the divergence theorem we have that $$\tag{1}\int_{\partial\Omega} (u\cdot\nu)v_{|\partial\Omega}=\int_\Omega v\operatorname{div}u+\int_\Omega \nabla v\cdot u $$
On the other hand, $$\tag{2}\langle \gamma (u),v_{|\partial\Omega}\rangle_{W^{-1+1/p,p},W^{1-1/p',p'}}=\int_\Omega \ell_p(v_{|\partial\Omega})\operatorname{div}u+\int_\Omega \nabla(\ell_p (v_{|\partial\Omega}))\cdot u $$
From $(1)$ and $(2)$, we conclude that (note that $v-\ell_p(v_{|\partial\Omega})\in W_0^{1,p'}(\Omega))$ $$\langle\gamma(u),v_{|\partial\Omega}\rangle=\int_{\partial\Omega} (u\cdot\nu)v_{|\partial\Omega},\ \forall\ v\in C^\infty(\overline{\Omega})$$
From the last equality, we conclude that $\gamma(u)=u\cdot\nu$ for every $w\in C^\infty(\overline{\Omega})^N$.