Gibbons-Hawking term and continuity The Gibbons Hawking boundary term is used to make a consistent variational procedure on manifolds with a boundary, and the usual form is given by:
$$\displaystyle \frac{\epsilon}{8\pi G}\int\limits_{\partial\mathcal{M}}d^3y\sqrt{|h|}K$$
where $K=n^\mu_{;\mu}$, and $\epsilon=\pm 1$ depending whether the surface is spacelike or timelike. Now, this all works nicely with spacelike/timelike surfaces, however if the surface happens to be lightlike (null), then the determinant $h=0$, but the metric tensor $g$ is also degenerate so the normal $n^\mu$ doesn't appear to be well defined, and there's no reason that I can see why the above boundary condition should also work in the null case. 
The most interesting case of a null surface is of course the event horizon of a black hole. The issue I have with this is that all the derivations of black hole entropy I've seen (including a PSE answer) just use the above term, with a slightly displaced boundary, and then in the end of the calculation simply take the zero limit of this displacement. My question is then, is this completely justified? Intuitively, since the event horizon is a coordinate singularity, and nothing "interesting" happens to an observer crossing it, everything should be nice and continuous, including the metric and its derivatives, so it seems plausible, but null surfaces seem to be a bit more involved mathematically. Basically, I can't come up with a convincing proof why the limiting procedure gives the exact same result as the one we'd obtain if we just started with a null surface from the beginning, I'm probably missing something obvious. Any help/references for this would be nice.
 A: I am not ready to give you an exact answer either, this is something I am working on myself (but I do so for the purpose of junction conditions, not BHs), however I can give you some pointers to consider.
The GHY term does not necessarily have to involve $K$. It comes from the fact that during the variation of the EH action, there is a term proportional to $g^{\mu\nu}\delta R_{\mu\nu}$, which when massaged further gives $\nabla_\mu Q^\mu$, where $Q$ is some vector field. We then use Gauss' theorem to turn it into a boundary integral. It can be shown that said boundary integral can expressed with $K$.
Gauss' theorem is $$ \int_D\text{div}X\sqrt{-g}d^nx=\int_ {\partial D}\langle n,X\rangle\sqrt{|h|}d^{n-1}\xi ,$$ where the angular brackets denote the metric tensor in an invariant sense and $h$ is the determinant of the induced metric and $\xi$ are hypersurface coordinates.
This theorem comes from Stokes' theorem the following way:
Stokes' theorem is $$ \int_D d\omega=\int_{\partial D}i^*(\omega), $$ where $\omega$ is an $n-1$ form and $i$ is the inclusion $i:\partial D\rightarrow M$. Let $\mu$ be the volume form associated with $g$, eg. one that is locally given by $\mu=\sqrt{-g}dx^1\wedge...\wedge dx^n$. The divergence can be defined as $\text{div}X\mu=\mathcal{L}_X\mu$.
You can check that if $\mu$ is the Riemannian volume element and $\nabla$ is the Levi-Civita connection, then $\text{div}X=\text{Tr}\nabla X$ (normally divergence with respect to volume form and divergence with respect to a linear connection are different things, but for Riemannian manifolds they agree).
With this we have $$ \int_D \text{div}X\mu=\int_D \mathcal{L}_X\mu,$$ but then we recall Cartan's formula as $$ \mathcal{L}_X=d\circ i_X + i_X\circ d, $$ where $i_X$ is the interior product (contraction with $X$). Since $\mu$ is closed we have $$ \mathcal{L}_X\mu=di_X\mu, $$ so by Stokes' theorem $$ \int_D \mathcal{L}_X\mu=\int_D d\ i_X\mu=\int_{\partial D}i^*(i_X\mu) $$ (execure me for using the same letter for the inclusion and the interiour product).
Now, one may show that if $X$ is decomposed as $X_n n+X^ie_{(i)}$, then $$ i_X\mu(Y_1,...,Y_{n-1})=\mu(X,Y_1,...,Y_{n-1})=X_n\mu(n,Y_1,...,Y_{n-1})=X_n i_n\mu(Y_1,...,Y_{n-1}), $$ because, due to the pullback $i^*$, I have restricted $i_X\mu$ to only act on $\partial D$-tangent vectors (so $Y_i$ are all tangent to the boundary), and the $e_{(i)}$ vector fields form some frame for the boundary, so because a differential form will only give nonzero results on linearly independent vectors, only the $X_n$-part of $X$ survives this plugin.
You can then show, that $i_n\mu$ is actually the same as $\sqrt{|h|}d^{n-1}\xi$.

The point is, instead of using the usual version of Gauss' theorem, you can use $\int_D\text{div}X\ \mu=\int_{\partial D}i^*(i_X\mu)$ instead, and you can try to choose a transverse vector field $N$, instead of a normal vector field $n$, and do a non-orthogonal decomposition of $X$ as $X=X_N N+X^ie_{(i)}$, and see what you get.
Maybe you can massage the expression so that it will involve some modification of $K_{\mu\nu}$. I would list these possible modifications but this post is getting too long. In the junction condition literature, these are known. Check out this article: https://arxiv.org/abs/gr-qc/0201054, here the geometry of null hypersurfaces and junction conditions along null hypersurfaces are treated in great detail. You may try to use some alternate quantity instead of $K$ along with what I detailed above to give a different form of the GHY boundary term, one that works for null surfaces as well.
