You should pull back the Riemannian metric $g$ on $X$ under the map $\alpha$. This defines a Riemannian metric $\alpha^*g$ on $\mathbb R^n$. The norm of $\nabla (f\circ \alpha)$ in this metric gives the answer to your question.
I'll illustrate this by an example. Let $X$ be the unit sphere in $\mathbb R^3$ with the metric induced from $\mathbb R^3$. The map $\alpha$ is $$\alpha(\theta,\phi)=(\cos\theta \sin \phi, \sin\theta\sin\phi, \cos\phi)\tag1$$
What does it mean to pull back the metric under $\alpha$? It means that the length of vector $p\frac{\partial}{\partial \theta}+q\frac{\partial}{\partial \phi}$ from the tangent space to $\mathbb R^2$ at a point $(\theta,\phi)$ is the length of
$$\left(p\frac{\partial}{\partial \theta}+q\frac{\partial}{\partial \phi}\right)\alpha = (-\sin\theta \sin\phi, \cos\theta\sin\phi,0)\,p+
(\cos\theta \cos \phi, \sin\theta\cos\phi, -\sin\phi)\,q
$$
in the metric of $X$ (which, for me, comes from $\mathbb R^3$). Using the Pythagorean theorem, I calculate
$$\bigg\|p\frac{\partial}{\partial \theta} + q\frac{\partial}{\partial \phi}\bigg\|_{\alpha^*g}^2 = (\sin^2\phi)\, p^2 +q^2 \tag2$$
Now suppose I want to calculate the norm of gradient of the function $f(x,y,z)=xyz$. Working through $\alpha$, I find
$$f\circ\alpha(\theta,\phi)=\cos\theta\sin\theta \sin^2\phi\cos\phi $$
and then calculate the gradient of $f\circ \alpha$ from its partial derivatives
$$\nabla(f\circ \alpha) = (\cos 2\theta \sin^2\phi\cos\phi)\frac{\partial}{\partial \theta} + \cos\theta\sin\theta (2\sin\phi\cos^2\phi-\sin^3\phi)\frac{\partial}{\partial \phi} $$
According to (2), the norm of the gradient is the square root of
$$ \cos^2 2\theta \sin^6\phi\cos^2\phi + \cos^2\theta\sin^2\theta (2\sin\phi\cos^2\phi-\sin^3\phi)^2 $$
Ugh, but this how it works in general.
If your parametrization $\alpha$ is conformal, the computations simplify because the pullback of $g$ under $\alpha$ is a multiple of the standard metric on $\mathbb R^n$. However, most $n$-manifolds with $n>2$ do not admit conformal parametrization.
First of all you have to note that
$$\hbox{div}(X)=\sum_{j=1}^{n}\langle\nabla_{E_j}X,E_j\rangle$$
taking in account that $E_i$ are orthonormal.
There fore by putting $X=\sum_{i=1}^{n}f_iE_i$, we can easily get
$$\hbox{div}(X)=\sum_{i,j=1}^{n}\langle\nabla_{E_j}(f_iE_i),E_j\rangle
=\sum_{i,j=1}^{n}\bigl(E_j(f_i)\langle E_i,E_j\rangle-f_i\langle\nabla_{E_j}E_i,E_j\rangle\bigr)=\sum_{i=1}^{n}E_i(f_i)$$
Best Answer
It's been a while, so sorry if this may be a bit too complicated:
One characterization of $\nabla f$ (which I use for the gradient of $f$) is that
$$g(Y,\nabla f)= df(Y) = Yf$$ for every $Y\in TM$.
If you assume now $Xf = |X|_g^2 = g(X,X)$, you get \begin{equation}\tag{*} \quad g(X, \nabla f) = g(X,X) \end{equation} Since $f$ is constant on it's level sets, for every curve $c:(-\varepsilon, \varepsilon) \rightarrow f^{-1}(c)$ the function $g(t) = f(c(t))$ is constant. Every tangent vector $Y = c^\prime(0)$ to the level set $f^{-1}(c) $ can be written that way. Then
$$0 = g^\prime(0) = df(c(0))c^\prime(0)= Yf|_{c(0)} = g(Y, \nabla f)|_{c(0)}$$ Since $Y$ is an arbitrary tangent vector to the level set, this implies that the gradient of $f$ is orthogonal to the level sets. Since you know by assumption that $X$ is also normal to the level sets, $\nabla f$ is a multiple of $X$, i.e. $\nabla f = \lambda X$. Now with $(*)$ you get $\lambda \equiv 1$