In this paper the following relation can be found:
\begin{align} \nabla\phi=\mathbf{g}^j\frac{\partial \phi}{\partial \xi^j}= J\frac{\partial}{\partial \xi^j}\left(\frac{\mathbf g^j}{J}\phi\right) \end{align}
I don't see how we went from the first form to the second. I've tried proving it backwards, i.e. I used the product rule on the second form and this led to: $$\mathbf{g}^j\frac{\partial \phi}{\partial \xi^j}+\phi\left(\frac{\partial \mathbf{g}^j}{\partial \xi^j}+\mathbf{g}^j\frac{\partial}{\partial \xi^j}\left(\frac 1J\right)\right)$$
The two terms in parentheses should cancel each other out. It vaguely resembles Ricci's theorem: $$\mathbf g^i_{,j}+\mathbf g^l\Gamma^i_{lj}=\mathbf 0$$
But I am completely stuck at this point. How can I proceed? Is there some other straightforward way of proving the first relation?
P.S.: I am aware there are other ways to represent the transformation of vector operators (using Lamé coeffecients $h_i$ for instance). But what I want to know is how the Jacobian snuck into this expression.
Note that $$\tag{1} J\frac{\partial}{\partial \xi^j}\left(\frac{\mathbf g^j}{J}\phi\right) = \frac{\partial\phi}{\partial \xi^j}\mathbf g^j + \phi J \frac{\partial }{\partial \xi^j}\left(\frac{1}{J}\mathbf g^j\right), $$
where the Einstein summation convention is applied.
The expression on the LHS of (1) is known as the conservative form of the gradient in terms of contravariant basis vectors. It is equivalent to the form you are expecting, which is
$$\nabla \phi = \frac{\partial\phi}{\partial \xi^j}\mathbf g^j,$$
so it remains to prove that
$$\tag{2}\phi J \frac{\partial }{\partial \xi^j}\left(\frac{1}{J}\mathbf g^j\right)= 0.$$
In the notation of the linked paper, $J^{-1}$ is the Jacobian determinant $\displaystyle\left|\frac{\partial x^i}{\partial \xi^j}\right|$ where $(x^1,x^2,x^3)$ are Cartesian coordinates.
Denoting the position vector as $\mathbf{r} = x^i \mathbf{e}_i$ (again using the Einstein convention) in terms of Cartesian coordinates and basis vector, the covariant and contravariant basis vectors in the curvilinear coordinate system are
$$\mathbf{g}_i = \frac{\partial \mathbf{r}}{\partial \xi_i} = \frac{\partial x_k}{\partial \xi_i}\mathbf{e}_k, \quad \mathbf{g}^i = \nabla \xi^i = \frac{\partial \xi^i}{\partial x_k}\mathbf{e}_k$$
The contravariant basis vectors are related to the covariant basis vectors by
$$\tag{3}\mathbf{g}^i = J \mathbf{g}_j \times \mathbf{g}_k,$$
where $(i,j,k)$ is a cyclic permutation, i.e. $(1,2,3)$, $(2,3,1)$, or $(3,1,2)$.
We can write the summed partial derivatives in (2) as
$$\tag{4}\sum_{i=1}^3\frac{\partial}{\partial \xi^i}(J^{-1}\mathbf{g}^i) = \sum_{cyc}\frac{\partial}{\partial \xi^i}(\mathbf{g}_j \times \mathbf{g}_k)\\ = \sum_{cyc}\left(\frac{\partial\mathbf{g}_j}{\partial \xi^i} \times \mathbf{g}_k +\mathbf{g}_j\times \frac{\partial\mathbf{g}_k}{\partial \xi^i}\right)$$
It is now easy to show that the RHS of (4) vanishes since we get canceling pairs of terms like
$$\frac{\partial\mathbf{g}_2}{\partial \xi^1} \times \mathbf{g}_3 = \frac{\partial^2\mathbf{r}}{\partial \xi^1 \partial \xi^2} \times \mathbf{g}_3,$$ and
$$\mathbf{g}_3 \times \frac{\partial\mathbf{g}_1}{\partial \xi^2} = \mathbf{g}_3 \times \frac{\partial^2\mathbf{r}}{\partial \xi^2 \partial \xi^1} = - \frac{\partial\mathbf{g}_2}{\partial \xi^1} \times \mathbf{g}_3$$
Thus, equation (2) holds.