Let $(M,G)$ be a Riemannian manifold of dimension $n$ and $S$ be a submanifold of $M$ of dimension $m$. Consider a vector field $X$ on $S$. $X$ is also a vector fieldd on $M$, so we can compute its divergence in some coordinates chart $x$ for $M$ as $$\text{div}X = \frac{1}{\sqrt{D}}\sum_{i = 1}^n \frac{\partial}{\partial{x}_i} \left( \sqrt{D} X^i \right)$$ with $D = |\det{G}|$.
Let now $g$ be the metric on $S$ induced by the one on $M$ (that is, $g$ is the pull-back of $G$ along the inclusion).
$X$ is a vector field on $S$, so we can compute its divergence in some coordinates chart $y$ for $S$ as $$\text{div}X = \frac{1}{\sqrt{d}}\sum_{i = 1}^m \frac{\partial}{\partial{y}_i} \left( \sqrt{d} X^i \right)$$ with $d = |\det{g}|$.
Is there any relation between the two expressions? In particular, if $X$ is divergence-free as a vector field in the ambient manifold, is it also divergence-free as a vector field on the submanifold?
The answer seems to be "yes" from an example, but I'd like to get a more theoretical understanding.
Example Let $M$ be the three-dimensional Euclidean space and $S$ the standard sphere. The vector field $$(y+z) \partial_x + (-x+z)\partial_y + (-x-y)\partial_z$$ is parallel to the sphere and divergence-free with respect to the Euclidean metric, and after a calculation it's divergence with respect to the pull-back metric is $$\sqrt{2} \left( \frac{1}{\tan{(\theta)}} - \frac{1}{\tan{\left(\theta \right)}}\right) \sin{\left(\phi + \frac{\pi}{4} \right)} = 0 $$