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.
