I used Gauss' identity to derive
$$(\ast) {1\over \varepsilon^{n-1}}\int_{\partial B(a,\varepsilon)} f dS = {1\over r^{n-1}}\int_{\partial B(a,r)} f dS $$
where $0<\varepsilon<r$ and $f$ is harmonic and $B(a,r)$ is the ball of radius $r$ around $a$.
Now I want to use this to prove that
$$ f(a) = {1\over |\partial B(a,r)|}\int_{\partial B(a,r)} f dS$$
Hence my question:
How to calculate $\lim_{\varepsilon \to 0}{1\over \varepsilon^{n-1}}\int_{\partial B(a,\varepsilon)} f dS$?
Does $ \lim_{\varepsilon \to 0}{1\over \varepsilon^{n-1}}\int_{\partial B(a,\varepsilon)} f dS = f(a)$? (I'm not sure what the value of $|\partial B(a,r)|$ is in higher dimensions, I believe it is the surface are of the sphere so I am also missing a factor on the RHS)
When scaled by the surface area, you do indeed get $f(a)$. Roughly speaking this is because you have an average of values all of which are being made arbitrarily close to $f(a)$. Precisely speaking you can say essentially the same thing: given $\varepsilon$, find $\delta$ so that if $|x-a|<\delta$ then $|f(x)-f(a)|<\varepsilon$. Let $g(x)=f(x)-f(a)$. Then for $r<\delta$, the average of $f$ over the sphere is $f(a)$ plus the average of $g(x)$, which is at most $\varepsilon$ in magnitude. So the average of $f$ is within $\varepsilon$ of $f(a)$.
Scaling only by $r^{n-1}$ as you are, let $\sigma_n$ be the surface area of the unit sphere. Then I'll show the limit is $\sigma_n f(a)$ as follows:
$$\left |-\sigma_n f(a)+\frac{1}{r^{n-1}} \int_{\partial B(a,r)} f(x) dS(x) \right | = \left | -\sigma_n f(a) + \frac{\sigma_n}{|\partial B(a,r)|} \int_{\partial B(a,r)} (f(a) + g(x)) dS(x) \right | \\ = \left | -\sigma_n f(a) + \frac{\sigma_n |\partial B(a,r)| f(a)}{|\partial B(a,r)|} + \frac{\sigma_n}{|\partial B(a,r)|} \int_{\partial B(a,r)} g(x) dS(x) \right | \\ \leq \frac{\sigma_n}{|\partial B(a,r)|} \int_{\partial B(a,r)} |g(x)| dS(x) \\ \leq \frac{\sigma_n \varepsilon |\partial B(a,r)|}{|\partial B(a,r)|}$$
provided $r<\delta$ as defined in the first paragraph.