First observe that on a compact Riemann manifold $(M, g)$ the operator $1+t \Delta$, $t>0$. $\Delta = d^*d: C^\infty(M)\to C^\infty(M)$ has a unique fundamental solution. Jacques Hadamard has constructed very explicit asymptotic expansions for this fundamental solution which lead to convergent series in the case of real analytic manifolds and metrics.
A modern description of Hadamard's construction can be found in volume 3, Sec. 17.4 of L. H"ormander's four volume on linear partial differential operators. (Hadamard's original memoir is also very useful, but harder to penetrate.) I will give a brief description of the fundamental solution $S_r$ of $(r+\Delta)$, $r>0$.
For $\nu=0,1,2,\dotsc $ and $r>0$ denote by $F_{\nu,r}(x)$ the generalized function (a.k.a. distribution) on $\mathbb{R}^n$ described as $\newcommand{\ii}{\boldsymbol{i}}$ a Fourier transform of a temperate distribution.
$$F_{\nu,r}(x)= \nu! (2\pi)^{-n} \int_{\mathbb{R}^n} e^{ \ii\langle x,\xi\rangle} (|\xi|^2+r)^{-\nu-1} d\xi. $$
The function $F_\nu$ can be expressed explicitly in terms of Bessel functions. Note that
$$F_{\nu,r}(x)= r^{\frac{n-\nu-1}{2}} F_\nu(\sqrt{r} x),\;\;F_\nu(x):=F_{\nu,r=1}(x). $$
If $$\Delta=-\sum_k\partial^2_{x_k}$$
denotes the (geometers') Laplacian in $\mathbb{R}^n$$ then
$$(r+\Delta)F_{0,r}=\delta_0,\;\;(r+\Delta)F_{\nu,r}=\nu F_{\nu-1, r},\;\;\forall \nu>0. $$
One can show that the generalized function $F_\nu$ depends only on the distance $|x|$.
Going back to the Riemann manifold $(M,g)$ we denote by $d: M\times M\to \mathbb{R}$ the geodesic distance function.
The Green function $G(x,y)$ then has an asymptotic expansion
$$ G(x,y)\sim \sum_{\nu=0}^\infty U_\nu(x,y) F_{\nu,r}( \; d(x,y)\;)$$
valid for $d(x,y)$ sufficiently small, where the functions $U_\nu(x,y)$ are explicitly described in the above reference. If $(M,g)$ is real analytic, then the above series converges in an appropriate sense.
This asymptotic expansion ought to be enough to investigate your question.