$\newcommand{\dd}{\text{d}}$
In the proof of symmetry of Green's function, Evans uses this estimate. Green's function is defined as $G(x,y) := \Phi(y-x) + \phi^x$, where $\Phi$ is the fundamental solution for Laplace's equation, and $\phi^x$ is the corrector function, depending on the geometry and satisfying $\Delta \phi^x = 0$ in $\Omega$ and $\phi^x(y) = G(x,y)$ on $\partial \Omega$.
How does one obtain the factor $\varepsilon^{n-1}$ in the following estimate $$\bigg|\int_{ \partial B_\varepsilon(x) } \frac{\partial G(x,z)}{\partial \nu} \, G(y,z) \, \dd \sigma(z) \bigg| \leq C \, \varepsilon^{n-1} \, \lVert G(x,y) \rVert_{L^\infty(\partial B_\varepsilon(x))} \overset{ \varepsilon \to 0}{\to} 0$$
I remember from previous sections, that $|\nabla \Phi(x)| \leq C \, |x|^{1-n}$, but I am not sure if this helps.
Thanks for any thoughts on this!
I believe the source of your confusion is that you have mixed up the $x$ and $y$ in Evans' proof. The estimate he claims is $$ \left\vert \int_{\partial B(x,\epsilon)} \frac{\partial w}{\partial \nu} v dS\right\vert \le C \epsilon^{n-1} \sup_{\partial B(x,\epsilon)} |v| = o(1) $$ where (and here is the key part) $$ w(z) = G(y,z) \text{ and } v(z) = G(x,z). $$
Consequently, near $x$, the function $w(z) = G(y,z)$ is smooth, and so the first order derivatives are bounded. Thus we may assume that $\epsilon < \epsilon_0$ and bound $$ \left\vert \int_{\partial B(x,\epsilon)} \frac{\partial w}{\partial \nu} v dS\right\vert \le n \alpha(n) \epsilon^{n-1} \sup_{\partial B(x,\epsilon_0)} |\nabla w| \sup_{\partial B(x,\epsilon)} |v| \le C \epsilon^{n-1} \sup_{\partial B(x,\epsilon)} |v| = o(1). $$ Here the appearance of $\epsilon^{n-1}$ is simply due to the surface area of $\partial B(x,\epsilon)$. The fact that the resulting estimate is $o(1)$ relies on bounds for the fundamental solution to the Laplacian (as in the previous sections in Evans) and the fact that the corrector function is harmonic and hence smooth.