It is actually a comment and question to the answer of Robert McNees in the following post:
Explicit Variation of Gibbons-Hawking-York Boundary Term
In deriving the variation of the extrinsic curvature, I obtain an additional contribution associate to what he calls $c^{\mu}$. I could not figure out why it should vanish in the first place (hence before assuming any boundary conditions as Dirichlet).
Thus for me the variation of $K$ is given by (where $D_{\mu}c^{\mu}$ is a total derivative w.r.t. $\int_{\partial M} \sqrt{h}$): $$\delta K=−12K_{\mu\nu}\delta g_{\mu\nu}−\tfrac{1}{2}n_{\mu}(\nabla_{\nu}\delta g_{\mu\nu}−g_{\nu\lambda}\nabla_{\mu}\delta g_{\nu\lambda})+D_{\mu}c^{\mu} + n^{\nu} (\partial_{\mu}n_{\nu}) c^{\mu}$$
The last contribution stems from the fact that $H^{\mu\nu}H_{\nu\rho}=\delta^{\mu}_{\rho}-n^{\mu}n_{\rho}$ is degenerated and yields not simply $\delta^{\mu}_{\rho}$, hence the usual Christoffel symbol structure has some further terms containing the normal directions.
Basically, what I was doing is the following: $$ \int \sqrt{h} \nabla_{\rho} c^{\rho}= -\int \sqrt{h} \tfrac{1}{2}(h^{\mu\nu}-g^{\mu\nu})(\partial_{\rho}g_{\mu\nu})c^{\rho}=\int \sqrt{h} \tfrac{1}{2}n^{\mu}n^{\nu}(\partial_{\rho}g_{\mu\nu})c^{\rho}$$ and then playing around with the orthogonality of $h$ and $n$. In my case it was $\nabla_{\rho}c^{\rho}$ appearing in the variation of $K$ in the first place and I converted in the above way to some total derivative and the additional term.
I am mentioning this point (though this question is old and it might have a trivial solution) because this post appears as the first hit whenever one searches for Gibbons-Hawking variation. Thus, any clarification why this additional term should disappear except for Dirichlet conditions might be helpful, or in the unlikely case I truly haven't overlooked anything, I think this term should be added in the previous results. Anyway, I am open for suggestions to match these outcomes.
