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.
The Gibbons-Hawking boundary term for a spacetime manifold is explicitly,
$$S_{GH}=\frac{1}{8\pi G}\int_{\partial M} d^3x \, \sqrt{|h|} \, K$$
where $\partial M$ is the boundary of $M$, $K$ the extrinsic curvature, and $h$ the determinant of the metric on the boundary. Let us Wick rotate the Schwarzschild metric to,
$$ds^2 = \left( 1-\frac{2GM}{r}\right)d\tau^2 + \left( 1-\frac{2GM}{r}\right)^{-1} dr^2 + r^2 d\Omega^2$$
We must impose a radial cutoff $R > GM$. The normal vector on the boundary is given by,
$$n=-\sqrt{1-\frac{2GM}{r}}\frac{\partial}{\partial r}$$
with a minus sign since we require the outward pointing normal, which points into the bulk. The metric on the boundary is then given by,
$$ds^2=\left( 1-\frac{2GM}{R}\right)d\tau^2 + R^2d\Omega^2$$
The extrinsic curvature is simply the divergence of the normal:
$$K=\nabla_a n^a = \frac{1}{r^2}\partial_r (r^2 n^r) \biggr\rvert_{r=R}= -\frac{2}{R}\sqrt{1-\frac{2GM}{R}} - \frac{GM}{R^2} \frac{1}{\sqrt{1-\frac{2GM}{R}}}$$
Can you take the calculation from here?
Best Answer
Regarding the last question, the Palatini formulation doesn't require the boundary term whereas standard GR does because the two formulations are not strictly equivalent. If the action is just the Einstein-Hilbert action, they give the same field equations, true, but for more generic actions (e.g., involving combinations and powers of the Ricci scalar $\mathcal R$), they can give different field equations.
Since we expect general relativity to be an effective theory with these extra terms being generated by quantum effects, it should be possible in principle to distinguish between default GR or Palatini formulation. In practice we don't have the technological means to do so right now.