Consider a (closed) Riemann surface and let $G(x,y)$ be the Green function of the Laplace-Beltrami operator. We can informally identify $G$ with the two-point correlation function for the Gaussian random field:

$$G(x,y)=\left<\phi(x)\phi(y)\right>=\frac{1}{Z} \int \mathcal{D}\phi\;\phi(x)\phi(y) \exp\left(-\frac{1}{2}\int \left| \nabla\phi(z)\right|^2 dV_g(z)\right).$$ This is easy to vary with respect to the metric $g_{\mu\nu}$:

$$\delta G(x,y)=-\frac{1}{2}\int dV_g(z)\;\delta g^{ij}(z)\left<\phi(x)\phi(y) \nabla_i\phi(z)\nabla_j\phi(z)\right>_c,$$ where $\left<\right>_c$ is the connected correlation function (the disconnected diagram $x\leftrightarrow y,z\leftrightarrow z$ is cancelled by the variation of the partition function). We can use Wick's theorem to compute this, getting

$$\boxed{\frac{\delta G(x,y)}{\delta g^{ij}(z)}=-\nabla_{[i} G(x,z)\nabla_{j]} G(z,y)}$$

(derivatives w.r.t. $z$). This formula looks wickedly similar to the Hadamard variation formula for the variation of the boundary of the domain in flat space. Yet I haven't been able to find any mentions of this in the mathematical literature.

Furthermore, if we define the regularized Green's function (a.k.a. the Robin function) by $$G^R(x)=\lim_{y\to x} \left(G(x,y)-\frac{1}{2\pi}\ln d(x,y)\right),$$ where $d(x,y)$ is the local geodesic distance, then I'm conjecturing the following variational formula: $$\frac{\delta G^R(x)}{\delta g^{ij}(z)}=-\nabla_{i} G(x,z)\nabla_{j} G(x,z)-\frac{1}{4\pi}\nabla_i \nabla_j G(x,z).$$ The second term is motivated by the well-known formula for conformal variations (where it is $\frac{1}{4\pi}\delta_x(z)$) and seems to be necessary to cancel the second order pole in this variation. Edit: this guess turned out to be wrong, see my answer below.