Suppose I have some vector field $X$ in $\mathbb{R^3}$ and a submanifold $M$ such that $\dim(M)=2$. I want to calculate the flux of $X$ through $M$, $\Phi_X(M)$. I would have two ways to do it, either go by the definition of flux or use the divergence theorem: $$\Phi_X(M)=\int_Vdiv(X)dV$$ where $V$ is the solid enclosed by $M$. Now, suppose that $M$ is some kind of quadric surface so that the use of a parametrization, $\psi$, would the make the integration much easier (Polar, spherical, etc.). If $X$ is given in the cartesian coordinate system along with the euclidean metric then the divergence of $X$ is straightforward, but if I'm using a parametrization $\psi$ for $V$ then I would practically be using a different coordinate system because $\dim V=3$, and so in that case, the divergence of $X$ would not be so straightforward as can be seen in this question and answer:Divergence in curvilinear coordinates. However I know very little (nothing actually) of covariance, contravariance and tensors, and this answer is only strict to orthogonal coordinate systems, and I'm wondering what other options are there to compute the flux by the divergence theorem and not calculate the divergence in a separate coordinate systems.
For example, if I have $M=\{(x,y,z)\in \mathbb{R^3}: z=x^2+y^2, 0\leq z\leq 1\}$ and the vector field $X:\mathbb{R^3} \to \mathbb{R^{3\times 1}}, (x,y,z)\to \begin{bmatrix} xz \\ z \\ -\frac{z}{2}\end{bmatrix}$. Then the euclidean divergence of $X$ is $div (X)=z-1/2$. Using $\psi$ as a parametrization of $V$, such that $\psi=\begin{bmatrix} \psi^1 \\ \psi^2 \\ \psi^3\end{bmatrix}$, would the flux be:
$$\Phi_X(M)=\int_{\psi^{-1}(V)}\bigg(\psi^3-\frac{1}{2}\bigg)dV$$
such that $dV=\sqrt{det(G(\psi;(x,y,z)))}d\psi^1d\psi^2d\psi^3$, where $G(\psi;(x,y,z))$ is the gram matrix of $\psi$ at $(x,y,z)$?