Sure. By a smooth dyadic decomposition it suffices to show that convolutions of the form $$ \varepsilon^{-n} \int_\Omega \varphi\left(\frac{y-z}{\varepsilon}\right)\ dy$$ for $0 < \varepsilon \lesssim 1$ and $\varphi$ a fixed bump function are smooth on $\partial \Omega$ uniformly in $\varepsilon$ (multiply by $\varepsilon^2$ and sum over dyadic $\varepsilon>0$ for a suitably chosen $\varphi$ to recover the Newton potential). This is trivial for large $\varepsilon$, so we may assume $\varepsilon$ small.
The strategy here is to transform this expression to eliminate all negative powers of $\varepsilon$, as this is the only obstruction to non-uniformity.
Locally we may parameterize $\Omega$ as a half-space $\{ (y', y_n): y_n \geq f(y') \}$ for some smooth function $f$, and then for $z = (z',f(z'))$ and $\varepsilon$ small enough the above expression becomes $$ \varepsilon^{-n} \int_{{\bf R}^{n-1}} \int_0^\infty \varphi\left(\frac{(y'-z', f(y')-f(z')+t)}{\varepsilon}\right)\ dy' dt$$ which after a change of variables $y' = z'+\varepsilon w$, $t = \varepsilon s$ becomes $$ \int_{{\bf R}^{n-1}} \int_0^\infty \varphi\left(\left(w, \frac{f(z'+\varepsilon w)-f(z')}{\varepsilon}+s\right)\right)\ dw ds. \quad (1) $$ By the fundamental theorem of calculus we have $$ \frac{f(z'+\varepsilon w)-f(z')}{\varepsilon} = \int_0^1 w \cdot \nabla f(z' + \varepsilon \theta w)\ d\theta$$ which can then be seen to (locally) be a smooth function of $z'$ and $w$ uniformly in $\varepsilon$. From this it follows from repeated differentiation under the integral sign and the chain rule that the expression in (1) is a smooth function of $z'$ uniformly in $\varepsilon$, giving the claim.