In Carlos Kenig's 1986 paper Elliptic Boundary Value Problems on Lipschitz Domains, he uses (but does not prove) the following result which he says is "very simple".
Fix $n\geq 2$. Let $\lambda$ be the function satisfying $\lambda(0)=0$ and $\lambda'(t) = \frac{1}{(1+t^2)^{n/2}}$. Let $\varphi:\mathbb{R}^{n-1}\to \mathbb{R}$ be a Lipschitz function and let $f\in C^\infty_c(\mathbb{R}^{n-1})$. Then he claims that
\begin{align*} &\lim_{\epsilon\to 0 } \int_{\|x-y\|>\epsilon} \frac{\varphi(y)-\varphi(x)-\nabla\varphi(x)\cdot(y-x)}{(\|x-y\|^2+|\varphi(x)-\varphi(y)|^2)^{n/2}} f(x) \, dx\\ &= - \lim_{\epsilon\to 0 } \int_{\|x-y\|>\epsilon} \frac{1}{\|x-y\|^{n-1}}\lambda\left( \frac{\varphi(x)-\varphi(y)}{\|x-y\|} \right) \nabla f(x)\cdot(y-x) \, dx \qquad(*) \end{align*}
for each point $y\in \mathbb{R}^{n-1}$ where $\varphi$ is differentiable.
This looks scary, not simple really, but I gave it a shot before asking here.
Let $\epsilon>0$ and let $y\in \mathbb{R}^{n-1}$ be such that $\nabla\varphi(y)$ exists. Then \begin{align*} &\int_{\|x-y\|>\epsilon} \frac{\varphi(y)-\varphi(x)-\nabla\varphi(x)\cdot(y-x)}{(\|x-y\|^2+|\varphi(x)-\varphi(y)|^2)^{n/2}} f(x) \, dx \\ &= \int_{\|x-y\|>\epsilon} \frac{\varphi(y)-\varphi(x)-\nabla\varphi(x)\cdot(y-x)}{\|x-y\|^n\left(1+\left(\frac{\varphi(x)-\varphi(y)}{\|x-y\|}\right)^2\right)^{n/2}} f(x) \, dx \\ &= \int_{\|x-y\|>\epsilon} \frac{\varphi(y)-\varphi(x)-\nabla\varphi(x)\cdot(y-x)}{\|x-y\|}\frac{1}{\|x-y\|^{n-1}}\lambda'\left(\frac{\varphi(x)-\varphi(y)}{\|x-y\|}\right) f(x) \, dx. \end{align*}
With $y$ fixed, we define a function $g$ on $\mathbb{R}^{n-1}\setminus\{y\}$ by $$ g(x) = \frac{\varphi(x)-\varphi(y)}{\|x-y\|}. $$
Notice that for almost every $x\in \mathbb{R}^{n-1}$ we have $$ \partial _jg(x) = \frac{\partial_j\varphi(x)\|x-y\|-(\varphi(x)-\varphi(y))\frac{x_j-y_j}{\|x-y\|}}{\|x-y\|^2} = \frac{\partial_i\varphi(x)}{\|x-y\|} - (\varphi(x)-\varphi(y))\frac{x_j-y_j}{\|x-y\|^3} $$ which implies $$ \nabla g(x)\cdot(x-y) = \frac{\nabla\varphi(x)}{\|x-y\|}\cdot(x-y) - (\varphi(x)-\varphi(y))\frac{\|x-y\|^2}{\|x-y\|^3} = \frac{\nabla \varphi(x)(x-y) - \varphi(x) + \varphi(y)}{\|x-y\|}. $$
Hence $$ \nabla(\lambda\circ g)(x)\cdot(y-x) = \lambda'\left(\frac{\varphi(x)-\varphi(y)}{\|x-y\|}\right)\frac{\varphi(y)-\varphi(x)-\nabla\varphi(x)\cdot(y-x)}{\|x-y\|} $$ by Chain Rule.
Thus the left hand side of $(*)$ is $$ \lim_{\epsilon\to 0 } \int_{\|x-y\|>\epsilon} \frac{1}{\|x-y\|^{n-1}}\nabla(\lambda\circ g)(x)\cdot(y-x) f(x) \, dx $$ and the right hand side is $$ - \lim_{\epsilon\to 0 } \int_{\|x-y\|>\epsilon} \frac{1}{\|x-y\|^{n-1}}(\lambda\circ g)(x) \nabla f(x)\cdot(y-x) \, dx. $$
However, I am now stuck. I am guessing there should be an integration by parts to make the derivative on $\lambda\circ g$ disappear, get the negative sign and get a derivative on $f$. I can't make this work. Also I don't think I have used the fact that $\lambda(0)=0$ or that $\varphi$ is differentiable at $y$ yet.
