We can define the extrinsic curvature of a codimension-one surface as
$$K_{ab} = q_a^{\phantom{a}c} q_b^{\phantom{b}d} \nabla_c n_d,$$
where $n^d$ is the normal vector to the hypersurface and $q_{ab} = g_{ab} - n_a n_b$ is the induced metric (with $g_{ab}$ the ambient metric). Here and below I'll assume the ambient geometry is Riemannian so that everything is spacelike. Now, the above expression can also be written in terms of the Lie derivative as
$$K_{ab} = \frac{1}{2}\mathcal{L}_n q_{ab}.$$
For a surface of codimension $n$, the extrinsic curvature can be defined in terms of the induced metric as
$$K^c_{ab} = -q_a^{\phantom{a}d} q_b^{\phantom{b}e} \nabla_d q_e^{\phantom{d}c},$$
where now the induced metric can be written as
$$q_{ab} = g_{ab} - \sum_i X^{(i)}_a X^{(i)}_b,$$
where the sum is over $n$ mutually orthogonal unit vectors normal to the surface. It is relatively simple to show that for any such unit normal vector $X^a$,
$$X_c K^c_{ab} = q_a^{\phantom{a}c} q_b^{\phantom{b}d} \nabla_c X_d.$$
My question is the following: can the above contraction be rewritten in terms of the Lie derivative, as it can in the codimension-one case? Specifically, is it true that
$$X_c K^c_{ab} = \frac{1}{2} \mathcal{L}_X q_{ab}?$$
I've been trying to derive the above expression for a while, but I've had no success. I'd like to know if I've just been screwing up my algebra, or if it's just not true.
The formula you wrote in terms of the lie derivative can work, as long as the normal vectors are chosen such that their commutator remains normal to the surface. This is always possible, but is not the most general choice allowed for the normal vectors. The choice of how to extend the normal vectors off the surface affects this commutator, but the extrinsic curvature vector itself does not depend on this choice.
It was a little tricky to show it, but I think I have a way. I'll do it for two spacelike unit normals, $n^a$ and $m^a$, and it should be clear that this generalizes to higher codimension surfaces. The first piece of information you need is that the projection of $\nabla_a n_b$ is symmetric. This follows from the Frobenius theorem: the integrability condition for the system of one forms $n_a$ and $m_a$ is $$ \begin{aligned} n\wedge m\wedge dn&=0\\ n\wedge m\wedge dm&=0. \end{aligned} $$ Decomposing $dn\equiv 2\nabla_{[a}n_{b]}$ into parallel and perpendicular components as $dn = dn_\parallel+n\wedge\alpha+m\wedge\beta+\gamma n\wedge m, $ the above integrability condition gives $$ n\wedge m\wedge dn_\parallel = 0, $$ which is equivalent to the statement that $dn_\parallel$ vanishes. So we can freely use the fact that $q_a^d q_b^e\nabla_d n_e = q_{(a}^dq_{b)}^e \nabla_d n_e$. Next compute $$ \begin{aligned} -n_c \nabla_d q_e^c &= -n_c\nabla_d(\delta_e^c-n_en^c-m_em^c) \\ &= \nabla_dn_e+m_en_c\nabla_dm^c \end{aligned} $$ Here I've used that the covariant derivative of $\delta^a_b$ is zero, and since $n^a$ is unit normalized, $n_c\nabla_a n^c = \frac12\nabla_a n^2 = 0$. Since we are eventually going to project this answer on the $e$ and $d$ indices, the second term will go away, so we can focus on the first. Now we compute the projections, which is a bit messy: $$ \begin{aligned} q_a^dq_b^e\nabla_dn_e&=(\delta_a^d-n_an^d-m_am^d)(\delta_d^e-n_dn^e-m_dm^e)\nabla_d n_e\\ &=(\delta_a^d-n_an^d-m_am^d)(\nabla_dn_b-m_bm^e\nabla_dn_e) \\ &=\nabla_an_b-n_aa_b-m_bm^e\nabla_an_e-m_am^d\nabla_dn_b+m_bn_am^en^d\nabla_dn_e+m_am_bm^dm^e\nabla_dn_e \end{aligned} $$ where $a_a = n^d\nabla_d n_a= £_n n_a$ is the acceleration, and the second expression in terms of the Lie derivative is true since $n$ is normalized. We should keep in mind that this expression is actually symmetric in $a$ and $b$, but I won't write in the symmetrizers to avoid clutter. Due to this symmetry, the first three terms can be written as $$ \begin{aligned} \nabla_a n_b-n_a a_b-m_bm^e\nabla_an_e &=\frac12£_ng_{ab}-n_b£_nn_a-m_b£_nm_a + m_b n^e\nabla_e m_a \\ &= \frac12 £_nq_{ab} + m_b n^e\nabla_e m_a \end{aligned} $$ This is almost what we want. Adding the extra unwanted term here to the remaining three terms in the above expression, and using the symmetry on $a$ and $b$, these terms can simplify to $$ \begin{aligned} &\,m_b(n^e\nabla_e m_a - m^e\nabla_e n_a+ n_am^en^d\nabla_dn_e+m_am^em^d\nabla_dn_e) \\ =&\, m_b([n,m]_a -n_a n_en^d\nabla_dm^e +m_am_em^d\nabla_dn^e) \\ =&\, m_b([n,m]_a-n_an^e[n,m]_e-m_am^e[n,m]_e) \end{aligned} $$ where the commutators are $[n,m]^a = n^e\nabla_e m^a-m^e\nabla_e m^a$. The above steps make frequent use of $n^e\nabla_d n_e = m^e\nabla_d m_e = 0$ and $n^e\nabla_d m_e = -m_e\nabla_d n^e$.
Now, the term in parentheses is simply the parallel component of the commutator $[n,m]$. This parallel component can always be chosen to be zero; if it is, it means that the normal vectors extend off the surface to be tangent to $2$ dimensional submanifolds. However, the most generic choice for extending these normal vectors off the surface will in general have a parallel component to this commutator, leading to the more general formula
$$ n_c K^c_{ab} = \frac12 £_n q_{ab} + m_{(b} q_{a)c}[n,m]^c $$
This then allows for an alternative expression for the extrinsic curvature tensor, $$ K^c_{ab} = \frac12n^c£_n q_{ab}+n^cm_{(a}q_{b)d}[n,m]^d+ \frac12 m^c£_mq_{ab}+m^cn_{(a}q_{b)d}[m,n]^d $$ and of course the generalization to higher codimension should be obvious.