I'm afraid that what I'm about to say is not quite you what to hear: your Neumann PDE does not have a solution for arbitrary choices of $g$.
To see why, let's compute the integral of $g$ over the surface $\Gamma$:
$$ \int_\Gamma g \ dS = \int_\Gamma \frac{\partial u}{\partial n} \ dS=\int_\Omega \nabla^2 u \ dV =\int_\Omega 0\ dV = 0.$$
Therefore, the average value of $g$ on the boundary surface $S$ must be zero. If $g$ fails to satisfy this condition, your PDE cannot be solved.
For a similar reason, your definition of the Green's function must be modified. Indeed, if we integrate $\frac{\partial G}{\partial n}$ over the surface $\Gamma$, we find that
$$ \int_\Gamma \frac{\partial G(\vec r, \vec r')}{\partial n_{\vec r}}dS(\vec r) = \int_\Omega \nabla_{\vec r}^2 G(\vec r,\vec r') \ dV(\vec r) = \int_V \delta(\vec r-\vec r') \ dV(\vec r) = 1.$$
So it is inconsistent to set $\frac{\partial G}{\partial n}$ equal to zero on the boundary $\Gamma$.
Fortunately, all is not lost! We can redefine the Green's function $G$ so that it satisfies
$$ \begin{cases} \nabla_{\vec r}^2 G(\vec r,\vec r') = \delta(\vec r - \vec r') &{\rm \ \ on \ \ } \Omega ,
\\
\frac{\partial G(\vec r, \vec r')}{\partial n_{\vec r}}=\frac 1 A & {\rm \ \ on \ \ } \Gamma, \end{cases}
$$
where $A = \int_\Gamma dS$ is the area of the boundary surface $\Gamma$.
Now, Green's identity states that
\begin{multline}\int_\Omega \left( u(\vec r) \nabla_{\vec r}^2 G(\vec r , \vec r') - G(\vec r,\vec r') \nabla_{\vec r} u(\vec r) \right) \ dV(\vec r) \\ = \int_\Gamma \left( u(\vec r) \frac{\partial G(\vec r,\vec r')}{\partial n_{\vec r}} - G(\vec r,\vec r') \frac{\partial u(\vec r)}{\partial n_{\vec r}}\right) dS(\vec r).\end{multline}
Plugging in my new definition of $G(\vec r, \vec r')$, we see that any $u(\vec r)$ satisfying your PDE must satisfy
we can immediately see that any solution to the PDE must satisfy
$$ u(\vec r') = - \int_\Gamma G(\vec r,\vec r')g(\vec r) \ dS(\vec r) + c, \ \ \ \ \ \ (\ast )$$
where the constant $c$ is equal to $\frac 1 A \int_\Gamma u(\vec r) dS(\vec r)$, the average value of $u$ on $\Gamma$.
Having redefined the Green's function, I'll give you an explicit expression in the case where $\Omega$ is a two-dimensional circular disk of radius $1$. Here it is:
$$ G(\vec r, \vec r') = \tfrac 1 {2\pi} \ln | \vec r -\vec r' | + \tfrac 1 {2\pi} \ln |\vec r - \vec r''|,$$
where $\vec r'' = \vec r'/|\vec r'|^2$ is the image of $r'$ under an inversion about the unit circle. It is similar to the Dirichlet Green's function, expect that we have a plus sign in front of the "image" term instead of a minus sign.
I believe that if you plug this Green's function into $(\ast )$ (with the constant $c$ chosen arbitrarily), you do indeed get a solution to your PDE, provided that $g$ obeys the consistency condition $\int_\Gamma g \ dS = 0$. This is stated in these lecture notes from my university, and also in Riley, Hobson and Bence, though I would love to see a rigorous proof! I wonder if anyone knows a good reference?
Best Answer
Edit: I have changed the proof of the \eqref{cc}$ \implies $\eqref{np} to correct an error in the reasoning pointed out by the asker. The development is now necessarily more complicated but entirely correct. I'd like to thank Prof. Alberto Cialdea for the useful discussion on the topic and the suggestion to use Fredholm theory and the equivalent Neumann problem for the Laplace equation.
What we want to prove is that the following Neumann problem $$ \color{green}{ \begin{cases} \Delta \varphi(x)=f(x) & x\in\bar{\Omega}\\ \frac{\partial}{\partial \vec{n}}\varphi(x)=g(x)& x\in\partial\bar{\Omega} \end{cases}\label{np}\tag{NP}} $$ is solvable if and only if the following compatibility condition $$ \color{blue}{ \int\limits_\bar{\Omega}f(x)\mathrm{d}x=\int\limits_{\partial\bar{\Omega}}g(x)\mathrm{d}\sigma_x. \label{cc}\tag{CC}} $$ holds (with obvious meaning of the symbols), i.e. \eqref{np}$ \iff $\eqref{cc}. Let's proceed with proving the two opposite implications.
Final notes on the method of proof of the implication \eqref{cc}$ \implies $\eqref{np}.
[1] V. P. Mikhailov (1978), Partial differential equations, Translated from the Russian by P.C. Sinha. Revised from the 1976 Russian ed., Moscow: Mir Publishers, p. 396 MR0601389, Zbl 0388.3500.
[2] V. S. Vladimirov (1971)[1967], Equations of mathematical physics, Translated from the Russian original (1967) by Audrey Littlewood. Edited by Alan Jeffrey, (English), Pure and Applied Mathematics, Vol. 3, New York: Marcel Dekker, Inc., pp. vi+418, MR0268497, Zbl 0207.09101.