In a curved $(3+1)$ dimensional spacetime with metric components $g_{\mu \nu}$, the covariant derivative of a $4$ vector $\mathbf V = (V^0, \vec V)$ is given by $$\nabla_\mu~ V^\mu = \frac{1}{\sqrt{-g}}\partial_\mu (\sqrt{-g}~V^\mu).$$
I expect that this relation can also be used to derive the expression for divergence of the $3$ vector $\vec V$ in a flat spatial hypersurface in a curvilinear coordinate system, eg. the cylindrical polar coordinates $(r,\phi,z)$. We will then need to replace the $\sqrt{-g}~$ by $\sqrt g~$ as the metric of the spatial hypersurface has a positive determinant. This will then give $$\vec \nabla \cdot \vec V = \nabla_i V^i = \frac{1}{r}\partial_r(r~V^r) + \partial_\phi V^\phi + \partial_zV^z.$$
However, the actual expression for the divergence of a $3$ vector in cylindrical polar coordinates is $$\vec \nabla \cdot \vec V = \frac{1}{r}\partial_r(r~V^r) + \frac{1}{r}\partial_\phi V^\phi + \partial_zV^z.$$ Can you please point out and explain where I am going wrong?