This should follow from local boundary $H^2(\Omega)$ regularity, applied separately near $\Gamma_0$ and $\Gamma_1$. Since the two boundary conditions are prescribed on separated regions (which is what I assume you mean by an annular domain), you can localise to one and ignore the other.
In what follows, note that a weak solution $u$ satisfies the weak formulation
$$ \int_{\Omega} \nabla u \cdot \nabla \varphi + u \varphi \,\mathrm{d}x = \int_{\Omega} f \varphi \,\mathrm{d}x + \int_{\Gamma_1} (g-u)\varphi \,\mathrm{d}x$$
for all $\varphi \in H^1_{\Gamma_0}(\Omega)$, that is all $\varphi \in H^1(\Omega)$ vanishing on $\Gamma_0$.
For regularity near $\Gamma_0$, we can take test functions $\varphi \in H^1_0(\Omega)$ and observe that the $\Gamma_1$ term vanishes in the weak formulation. Then $H^2$ regularity follows in the standard way by:
Locally flattening the boundary near some $x_0 \in \Gamma_0$,
Taking tangential difference quotients to show all the tangential derivatives $\partial_iu$ locally lie in $H^1$ for $1 \leq i \leq n-1$. Here it is convenient that the Dirichlet boundary is zero, which ensures the difference quotients remain valid test functions.
Using the equation deduce that the final derivative $\partial_nu$ also locally lies $H^1$.
This is standard, and can be found in any graduate level PDE text covering elliptic equations. See for instance Chapter 6, Section 3.2 of
Evans, Lawrence C., Partial differential equations, Graduate Studies in Mathematics. 19. Providence, RI: American Mathematical Society (AMS). xvii, 662 p. (1998). ZBL0902.35002.
For regularity near $\Gamma_1$, the strategy is similar, involving flattening the boundary and using difference quotients. Here $\varphi$ built from difference quotients will automatically be a valid test function (assuming your cutoff vanishes near $\Gamma_0$), but you need to estimate an additional term coming from the $\Gamma_1$ integral. I will leave this an exercise, but since you are only taking difference quotients in the tangential directions it's just a matter of using the trace theorem to estimate $u$.
Best Answer
As suggested in the comments, for a fixed $k \in \mathbb R$ you can always find a weak (or variational) solution of the problem $$ \left\lbrace \begin{array}{cccl} -\Delta u + k &=& f & \text{in}& \Omega\\ \partial_nu &=& 0 & \text{on }& \partial \Omega \end{array}\right. $$ Note that you should impose regularity conditions on the boundary of $\Omega$ to that the trace is well defined, for instance Lipschitz boundary. It is clear that the solution is unique, hence you get a map $\Phi : \mathbb R \rightarrow L^2(\Omega)$.
Now you wish to use regularity theory in the sens that any $u \in H^1(\Omega)$ that solves the above problem is in fact $H^2$ with $$\|u\|_{H^2(\Omega)} \leq c \|\Delta u \|_{L^2(\Omega)} $$ where $c$ does not depend on $u$. To this aim you need to ask more regularity for $\partial \Omega$, I think that $\partial \Omega$ of class $C^2$ would be enough. If this is the case, $\Phi$ is in fact valued in $H^2(\Omega)$.
Finally we want to solve $$ \left\lbrace \begin{array}{cccl} -\Delta u + \int u &=& f & \text{in}& \Omega\\ \partial_nu &=& 0 & \text{on }& \partial \Omega \end{array}\right. $$ Note that this problem has at most one solution, so we are left to find a solution. Said in other words, we wish to find $u \in H^2(\Omega)$ such that $$ \Phi \int u = u. $$ Note that we can look at $u \in L^1(\Omega)$ solving the above equation: by regularity theory such a $u$ would be $H^2$. To this aim consider $\Psi : L^1(\Omega) \rightarrow L^1(\Omega)$ defined by $$ \Psi u = \Phi \int u, $$ it is enough to find a fixed point to $\Psi$. In view of the Schauder fixed point theorem and the Sobolev embeddings, it remains only to show that $\Psi$ is continuous. Being $\int : L^1(\Omega) \rightarrow \mathbb R$ continuous, it is enough to show that $\Phi : \mathbb R \rightarrow L^1(\Omega)$ is continuous. To this aim pick $\lambda,\mu \in \mathbb R$, write $\Phi(\lambda) = u$ and $\Phi(\mu) = v$ and compute $$ \|\Phi(\lambda) - \Phi(\mu)\|_{L^1(\Omega)} = \|u-v\|_{L^1(\Omega)} \lesssim \|u-v\|_{H^2(\Omega)} \leq c \| \Delta (u-v)\|_{L^2(\Omega)} = c |\lambda - \mu| |\Omega|^{1/2}. $$