In the case that you mentioned, we want to avoid this cut-off/difference quotients approach, since it could be hard to prove that $\partial_{x_i} (\xi \partial_{x_j} u)$ is a valid test function.
In general, when working with regularity theory, another standard approach is to use an 'approximated problem'. However, the kind of the approximated problem, of course, depends on the PDE.
For the Neumann like problem I suggest the following approximation:
First observe that since $\int_\Omega f = 0$, we can easily construct a sequence $\{f_n\} \subset C^\infty(\Omega)$ such that $f_n \to f$ in $L^2(\Omega)$ and $\int_\Omega f_n=0, \ \forall n \in \mathbb{N}$.
Then we consider $u_n \in C^\infty(\Omega)$ such that
$(*) -\Delta u_n + \frac{1}{n} u_n = f_n, \mbox{ in } \Omega \ \ $ and $\dfrac{\partial u_n}{\partial \nu}=0 \ , \mbox{on }\partial \Omega.$
The sequence $\{ u_n \}$ can be obtained by the use of Theorem 2.2.2.5, p.91 and Theorem 2.5.1.1, p. 121 of Grisvard's Book p. 91.
In fact you just need to use a boostrap argument to $-\Delta u_n = -\frac{1}{n} u_n + f_n$.
Notice that $\int_\Omega u_n=n\int_\Omega f_n = 0$.
Now, you use $u_n$ as your test functions and obtain the following estimate:
$(**) \|\nabla u_n\|_{L^2}^2 \leq \|f_n\|_{L^2}^2$, $\forall \ n \in \mathbb{N}$
Now you use $-\Delta u_n $, as a test function to your PDE ( observe that $-\Delta u_n$ is a valid test function, anyway we don't need to worry about it since the approximated equation holds everywhere).
After integrating by parts, by using $(**)$ with some standard manipulations with your boundary terms you end up with
$\|D^2 u_n \|_{L^2}^2 \leq C(\partial \Omega)\|f_n\|_{L^2}^2$, $\forall \ n \in \mathbb{N}$.
(For instance, see Grisvard's book p.132-138, in particular eq. 3.1.1.5)
The key point for the above estimation is to control the boundary elements in terms of the mean curvature of $\partial \Omega$.
Now, since $\int_\Omega u_n =0$ we conclude that
$\|u_n \|_{H^2}^2 \leq C(\partial \Omega)\|f_n \|_{L^2}^2 $
so that
$\|u_n \|_{H^2}^2 \leq C(\partial \Omega) \|f\|^2$, the $L^2$ norm of $f$.
In this way we obtain $u\in H^2$ such that $u_n \to u$ weakly in $H^2$ and strongly in $H^1$.
Observe that the latter convergence is sufficient to handle the term $\dfrac{1}{n}u_n$.
Then, we can pass to the limit in equation (*) so that $u$ is a strong solution of
$-\Delta u = f$ in $\Omega$
$\dfrac{\partial u}{\partial \nu}=0$ on $\partial \Omega$
with
$\|u\|_{H^2} \leq C(\partial \Omega) \|f\|$, the $L^2$ norm of $f$.
There exists a theory of Green functions for general parabolic boundary value problems which covers the case you are interested in, in particular papers by Eidelman, Ivasishen, Solonnikov. For references see
S. D. Eidelman and N. V. Zhitarashu, Parabolic boundary value problems. Basel: Birkhäuser (1998).
Unfortunately, most of the papers on this subject are available only in Russian.
Best Answer
Set $G$ a primitive of $g$. Then the solution is a critical point of the functional $E:H^1(\Omega)\rightarrow{\mathbb R}$ defined by $$E[u]:=\int_\Omega \left(\frac12|\nabla u|^2+fu\right)dx-\int_{\partial\Omega}G(u)ds.$$ If $G(\pm\infty)=-\infty$, you may look for a minimum of $E$ over $H^1(\Omega)$. When $g$ is a decreasing function, $G$ is concave and $E$ is strictly convex: your solution is unique.