The normal bundle along the submanifold in the Riemannian background is a special case of the following situation.
Let $M$ be a suitably good manifold of dimension $n$ with a Riemannian metric $g$. It has the Levi-Civita connection $\nabla = \nabla^g$ compatible with $g$, that is $\nabla g = 0$. Taking care with the orientation, there is also the Riemannian volume form $\mathrm{d}M$ on $M$.
Furthermore, let $E \to M$ be a (finite-dimesional) vector bundle over $M$ and forgive me if I sloppily denote the space of its sections by $E$ as well. Let $h: E \otimes E \to \mathbb{R}$ be a vector-bundle metric. Let $\nabla^E$ be some metric connection in $E$, that is $\nabla^E h = 0$.
There is a whole family of tensor products of the tangent, cotangent, and the given vector bundle $E$ on the manifold $E$, which are endowed with the corresponding tensor-product connection. More precisely, on each of such bundles, the connection acts by extending the Levi-Civita connection and the connection $\nabla^E$ using the Leibniz rule and linearity. Each of those connections we are going to denote simply by $\nabla$ and refer to them collectively as the coupled Levi-Civita connection.
The vector-bundle connection Laplacian $\Delta^E$ is defined on any section $X$ of the bundles, described above, by the formula
$$
\Delta^E X := g^{a b} \nabla_a \nabla_b X
$$
where $\nabla$ is the coupled Levi-Civita connection. The formula for the normal Laplacian in the OP is just a particular case of this.
Now I claim that
$$
\big( \nabla^a h(X, \nabla_a X) \big) \mathrm{d}M = \mathrm{d}\omega
$$
for some $(n - 1)$-form $\omega$ on $M$. This allows us to use the general Stokes theorem
$$
\int_{M} \mathrm{d}\omega = \int_{\partial M} \omega
$$
and, if we take into account the OP restrictions on the sections $X$ with regards to the boundary, we can say that
$$
\int_{M} \big( \nabla^a h(X, \nabla_a X) \big) \mathrm{d}M = 0
$$
Using the Leibniz rule, we calculate
$$
\nabla^a h(X, \nabla_a X) = h(\nabla^a X, \nabla_a X) + h(X, \nabla^a \nabla_a X) = |\nabla X|^2 + h(X, \Delta^E X)
$$
From the last display, integrating and using the Stokes theorem, we obtain
$$
\int_{M} |\nabla X|^2 \mathrm{d}M = - \int_M h(X, \Delta^E X) \mathrm{d}M
$$
Proof of the claim. This is going to be a calculation using an abstract index notation technique for dealing with differential forms.
A $(n)$-form $\xi$ using the index notation can be presented as $\xi_{a_1 \dots a_n}$, where the sequentially numbered indexes are assumed to be antisymmetrized:
$$
\xi_{a_1 \dots a_n} = \xi_{[a_1 \dots a_n]}
$$
The exterior derivative $\mathrm{d}$ on the differential form $\xi$ as above is defined by the formula:
$$
(\mathrm{d} \xi)_{a_0 a_1 \dots a_n} := (n + 1) \nabla_{a_0} \xi_{a_1 \dots a_n}
$$
Recall that the volume form $\mathrm{d} M$ is a ${n}$-form, so we can write it as $(\mathrm{d} M)_{a_1 \dots a_n}$ to emphasize that.
Now we want $\xi$ to be a $1$-form (a covector field). We can introduce a notation, just for the convenience sake, for the interior product of the covector $\xi$ with the volume form $\mathrm{d} M$ as so:
$$
(\iota_{\xi} \mathrm{d} M)_{a_2 \dots a_n} := g^{a a_1} \xi_a (\mathrm{d} M)_{a_1 \dots a_n}
$$
Clearly, $\iota_{\xi} \mathrm{d} M$ is a $(n-1)$-form.
Let us compute the exterior derivative of $\iota_{\xi} \mathrm{d} M$ to see what happens.
$$
\big( \mathrm{d} (\iota_{\xi} \mathrm{d} M) \big)_{a_0 a_2 \dots a_n} = n \nabla_{a_0 } (\iota_{\xi} \mathrm{d} M)_{a_2 \dots a_n}
= n \nabla_{a_0 } g^{a a_1} \xi_a (\mathrm{d} M)_{a_1 \dots a_n} = \\
= n g^{a a_1} \big( \nabla_{a_0 } \xi_a (\mathrm{d} M)_{a_1 \dots a_n} \big)
= n g^{a a_1} \big( (\nabla_{a_0 } \xi_a) (\mathrm{d} M)_{a_1 \dots a_n} + \xi_a \nabla_{a_0 } (\mathrm{d} M)_{a_1 \dots a_n} \big) = \\
= n g^{a a_1} \big( (\nabla_{a_0 } \xi_a) (\mathrm{d} M)_{a_1 \dots a_n}
$$
The last expression in the display above is a $(n)$-form on $M$, so it must be proportional to $\mathrm{d} M$:
$$
g^{a a_1} \big( (\nabla_{a_0 } \xi_a) (\mathrm{d} M)_{a_1 \dots a_n} = f (\mathrm{d} M)_{a_0 a_2 \dots a_n}
$$
Contracting with $\mathrm{d} M$ and using its properties (see R.Wald, General Relativity, p.433), we can find that
$$
\nabla^b \xi_b = n f
$$
We see that for any covector $\xi$ the so-called divergence term $(\nabla^a \xi_a) \mathrm{d} M$ is an exact differential form:
$$
(\nabla^a \xi_a) \mathrm{d} M = \mathrm{d}(\iota_{\xi} \mathrm{d} M)
$$
Taking $\xi_a = h(X, \nabla_a X)$ we obtain the claimed property.
