Explicit Variation of Gibbons-Hawking-York Boundary Term Are there any references that present the explicit variation of the Hilbert-Einstein action plus the Hawking-Gibbons-York boundary term, and demonstrate the cancellation of the normal derivatives of metric variations? I have tried to read the original papers by York and Gibbons&Hawking, but they are not that pedagogical to me.
 A: I am doing my Bachelors on something like that and we even chose to include double boundary terms in the variation.
In addition to this paper, Variational principle and 1-point functions in 3-dimensional flat space Einstein gravity by Stephane Detournay et al., we get one further term: 
$$-\frac{(1-\alpha)}{16\pi G}\int_{\mathcal{\partial M}}d^3x\sqrt{-\gamma}\frac{1}{2}n^c \nabla_c(n^an^b)\delta g_{ab} $$
the total variation of the GHY term (with a free parameter $\alpha$) and the ordinary EH term then looks:
\begin{align}
\delta \Gamma_{(\alpha)}&=\frac{1}{16\pi G}\int_{\mathcal{M}}d^4x\sqrt{-g}G^{ab}\delta g_{ab} \\
&\quad + \frac{1}{16\pi G}\int_{\mathcal{{\partial}M}}d^3x\sqrt{-\gamma}\left(K^{ab}-\alpha K g^{ab}+(2\alpha-1)Kn^an^b\right)\delta g_{ab}
\\&\quad+\frac{(1-\alpha)}{16\pi G}\int_{\mathcal{\partial M}}d^3x\sqrt{-\gamma}\left(\gamma^{ab}n^c\nabla_c-\frac{1}{2}n^c\nabla_c(n^an^b)\right)\delta g_{ab} \\&\quad+\frac{(2\alpha-1)}{16\pi G}\int_{\mathcal{\partial^2M}}d^2x\sqrt{-\gamma'}{n'}^an^b\delta g_{ab}
\end{align}
A: I've never seen a paper where the calculation is performed in a manifestly covariant manner. However, I've posted a set of reference notes on my website (http://jacobi.luc.edu/notes.html) that contains the variations needed to carry out the calculation. Let me summarize the calculation here.
The action for gravity on a compact region $M$ with boundary $\partial M$ is
$$I_{EH} + I_{GHY} = \frac{1}{2 \kappa^2} \int_{M}d^{d+1}x \sqrt{-g} R + \frac{1}{\kappa^2} \int_{\partial M} d^{d}x \sqrt{-h} K ~.$$
The metric on $M$ is $g_{\mu\nu}$, and $R = g^{\mu\nu} R_{\mu\nu}$ is the Ricci Scalar. The induced metric on the boundary $\partial M$ is $h_{\mu\nu} = g_{\mu\nu} - n_{\mu} n_{\nu}$, where $n^{\mu}$ is the (spacelike) unit vector normal to $\partial M \subset M$. Now consider a small variation in the metric: $g_{\mu\nu} \to g_{\mu\nu} + \delta g_{\mu\nu}$. The quantities appearing in the Einstein-Hilbert part of the action change in the following manner:
$$ \delta \sqrt{-g} = \frac{1}{2} \sqrt{-g} g^{\mu\nu} \delta g_{\mu\nu}$$
$$ \delta R = -R^{\mu\nu} \delta g_{\mu\nu} + \nabla^{\mu}\left(\nabla^{\nu} \delta g_{\mu\nu} - g^{\nu\lambda} \nabla_{\mu} \delta g_{\nu\lambda} \right)$$ 
Thus, the change in $I_{EH}$ is 
$$\begin{aligned}\delta I_{EH} = & \frac{1}{2\kappa^{2}}\int_{M} d^{d+1}x \sqrt{-g} \left(\frac{1}{2} g^{\mu\nu} R - R^{\mu\nu} \right)\delta g_{\mu\nu}\\
& + \frac{1}{\kappa^2} \int_{\partial M} d^{d}x \sqrt{-h} \frac{1}{2} n^{\mu} \left(\nabla^{\nu} \delta g_{\mu\nu} - g^{\nu\lambda} \nabla_{\mu} \delta g_{\nu\lambda}\right)~,\end{aligned}$$
with the boundary term coming from the volume integral of the total derivative in $\delta R$. The variations of the quantities in the GHY term are a bit more complicated to work out, but they all basically follow from standard definitions and this result for the variation of the normal vector:
$$\delta n_{\mu} = \frac{1}{2} n_{\mu} n^{\nu} n^{\lambda} \delta g_{\nu\lambda} = \frac{1}{2} \delta g_{\mu\nu} n^{\nu} + c_{\mu}~.$$
In the second equality I've introduced a vector $c_{\mu}$ that is orthogonal to $n^{\mu}$; it is given by
$$c_{\mu} = - \frac{1}{2} h_{\mu}{}^{\lambda} \delta g_{\nu\lambda} n^{\nu} ~.$$
The reason I've introduced this vector is that the variation in the trace of the extrinsic curvature can be written as
$$\delta K= - \frac{1}{2} K^{\mu\nu} \delta g_{\mu\nu} - \frac{1}{2} n^{\mu}\left(\nabla^{\nu} \delta g_{\mu\nu} - g^{\nu\lambda} \nabla_{\mu} \delta g_{\nu\lambda} \right) + D_{\mu} c^{\mu}$$
where $D_{\mu}$ is the covariant derivative along $\partial M$ that is compatible with the induced metric $h_{\mu\nu}$. So, the change in the GHY part of the action is
$$\delta I_{GHY} = \frac{1}{\kappa^2} \int_{\partial M} d^{d}x \sqrt{-h}\left(\frac{1}{2}h^{\mu\nu} \delta g_{\mu\nu} K + \delta K \right)~.$$ 
Combining this with $\delta I_{EH}$ we see that the several terms cancel, leaving
$$\begin{aligned} \delta I = & \frac{1}{2\kappa^2}\int_{M} d^{d+1}x \sqrt{-g}\left(\frac{1}{2} g^{\mu\nu} R - R^{\mu\nu} \right)\delta g_{\mu\nu}\\
& + \frac{1}{\kappa^2}\int_{\partial M} d^{d}x \sqrt{-h}\left(\frac{1}{2}(h^{\mu\nu} K - K^{\mu\nu})\delta g_{\mu\nu} + D_{\mu} c^{\mu} \right) ~.\end{aligned}$$
We discard the term $D_{\mu} c^{\mu}$, which is a total boundary derivative. 
