Say $(M,g)$ is a $2$-dimensional oriented Riemannian manifold, and let $\star$ be the Hodge star operator. Given a vector field $F$ on $M$, we define
\begin{align}
\text{curl}(F)&:=\star d(g^{\flat}F)
\end{align}
i.e we take our vector field $F$, convert it to a covector field $g^{\flat}(F)$, take its exterior derivative (so we now have a $2$-form) and finally take its Hodge dual to get a smooth function (since we're on a $2$-dimensional manifold). The fact that the curl of a vector field in $2$-dimensions yields a smooth function corresponds to your observation that there's only one non-vanishing term. The thing you're missing is the final Hodge star (the extra $r$ you have is the same $r$ in $dx\wedge dy=r\,dr\wedge d\theta$).
Explicitly, suppose we're in the plane and using polar coordinates. Suppose $F=f^r\mathbf{e}_r+ f^{\theta}\mathbf{e}_{\theta}$ is the expansion in terms of unit vector fields (as in vector calculus). Then, $g^{\flat}(F)=f^r\,dr +f^{\theta}\,rd\theta$ (here this factor of $r$ is the norm of $\frac{\partial}{\partial \theta}$, i.e $\sqrt{g_{\theta\theta}}=\sqrt{r^2}=r$). So,
\begin{align}
d(g^{\flat}(F)) &=
\left(\frac{\partial f^r}{\partial r}\,dr+\frac{\partial f^r}{\partial \theta}\,d\theta\right)\wedge dr +
\left(\frac{\partial (rf^{\theta})}{\partial r}\,dr+\frac{\partial (r f^{\theta})}{\partial \theta}\,d\theta\right) \wedge d\theta\\
&=\left(\frac{\partial (r f^{\theta})}{\partial r}-\frac{\partial f^r}{\partial \theta}\right)dr\wedge d\theta\tag{$*$}\\
&=\frac{1}{r}\left(\frac{\partial (r f^{\theta})}{\partial r}-\frac{\partial f^{r}}{\partial \theta}\right)\, rdr\wedge d\theta
\end{align}
In the last line, I simply divided and multiplied by $r$.
The reason for this is that the area-form in the plane is $dx\wedge dy=r\,dr\wedge d\theta$, and the nice thing about the Hodge star is that for any function $\phi$ and the area form $\omega$, we have $\star(\phi \, \omega)=\phi (\star\omega)=\phi\cdot 1$, i.e $\star(r\,dr\wedge d\theta)=1$. Therefore,
\begin{align}
\text{curl}(F)&=\star d(g^{\flat}(F))=\frac{1}{r}\left(\frac{\partial (r f^{\theta})}{\partial r}-\frac{\partial f^{r}}{\partial \theta}\right)
\end{align}
Where you went 'wrong'
Note that in your calculation, you stopped at $(*)$, i.e the expression you found is the component of a $2$-form (the coefficient of $dr\wedge d\theta$). Inherently, there's nothing wrong with that, but the classical definition of curl requires an extra Hodge star application. So really, it boils down to the fact that your initial expression $\nabla_av_b-\nabla_bv_a$ is the component of a $2$-form, but it is NOT the usual definition of curl.
Extra Ramblings.
For a vector field $F$ in $\Bbb{R}^3$, the definition of curl must be modified slightly:
\begin{align}
\text{curl}(F)&:=g^{\sharp}(\star d(g^{\flat}F))
\end{align}
i.e we convert the vector field $F$ to a $1$-form $g^{\flat}(F)$, take its exterior derivative (so we have a $2$-form now), then take the Hodge star (so we get a $3-2=1$-form) and finally use $g^{\sharp}$ to convert it back to a vector field (as a side remark: from this one can already see that the curl of a vector field is an EXTREMELY unnatural operation, and it's unfortunate we use it so often in elementary vector calculus).
BTW, if you want to think of a vector field in $\Bbb{R}^2$ as a vector field in $\Bbb{R}^3$, you can do so by thinking of it as having no $z$ component, in which case applying this 3-D definition of curl will yield a vector field having only a $\mathbf{e}_z$ component, which matches with what we found above (so overall there's no inconsistency in jumping between $2$ and $3$ dimensions).
Found the answer to my question as well as some new stuff. For starters, whilst $\widetilde{g}=\left(\nabla\vec{v}\right)\left(\nabla\vec{v}\right)^T$ is a correct equation, it isn’t what I’m looking for. The correct usage of this would be to describe the surface of a sphere in 3D with the coordinates u and v. The vector gradient works properly with this, and the resulting metric tensor is 2D and has the curvature of the sphere. You can see what’s actually happening here. We take a surface and merely change the coordinate system. But this cannot affect the surface’s curvature though. By default, we were remapping a flat 4D space to 4D space. So, this method cannot create a curved 4D metric with a 4D vector field. We could create any metric we wanted using a higher dimensional vector field due to the Nash embedding theorem, where the dimension required for an nD metric tensor is at most $\frac{n\left(3n+11\right)}{2}$, so 46D for a 4D metric tensor. Obviously, this is the antithesis of what I wanted.
The formula $\widetilde{g}=\nabla\vec{v}+\left(\nabla\vec{v}\right)^T$ has also been ruled out. The reasons why are as follows. The formula is linear, so
$$\widetilde{g}=\nabla\vec{v}+\left(\nabla\vec{v}\right)^T=\nabla\left({\vec{v}}_1+{\vec{v}}_2\right)+\left(\nabla\left({\vec{v}}_1+{\vec{v}}_2\right)\right)^T=\nabla{\vec{v}}_1+\nabla{\vec{v}}_2+\left(\nabla{\vec{v}}_1\right)^T+\left(\nabla{\vec{v}}_2\right)^T=\left(\nabla{\vec{v}}_1+\left({\vec{v}}_1\right)^T\right)+\left(\nabla{\vec{v}}_2+\left(\nabla{\vec{v}}_2\right)^T\right)={\widetilde{g}}_1+{\widetilde{g}}_2$$
All metric tensor fields should be possible, and their configurations determine the form of the stress energy tensor. As such, taking this and the linearity, we add various metric tensor fields together to get a metric tensor with only 1 non-zero unique term. We could create 10 of these fields and add them together to create any 4D metric tensor field. Now we need to check if this formula has the capacity to create such a field. We will start, and end, with the following case.
$$g_{11}=f\left(\vec{x}\right)=\sum_{m_1,m_2,m_3,m_4=0}^{\infty}{c_{m_1,m_2,m_3,m_4}x_1^{m_1}x_2^{m_2}x_3^{m_3}x_4^{m_4}}$$
$$v_1=\frac{1}{2}\int{f\left(\vec{x}\right)dx_1}$$
$$\left(\frac{\partial v_1}{\partial x_2}\right)=\frac{1}{2}\left(\frac{\partial}{\partial x_2}\int{f\left(\vec{x}\right)dx_1}\right)=\frac{1}{2}\int{\left(\frac{\partial f\left(\vec{x}\right)}{\partial x_2}\right)dx_1}=-\left(\frac{\partial v_2}{\partial x_1}\right)$$
$$v_2=-\frac{1}{2}\iint{\left(\frac{\partial f\left(\vec{x}\right)}{\partial x_2}\right){dx_1}^2}$$
$$0=\left(\frac{\partial v_2}{\partial x_2}\right)=-\frac{1}{2}\iint{\left(\frac{\partial^2f\left(\vec{x}\right)}{\partial x_2^2}\right){dx_1}^2}=-\frac{1}{2}\sum_{m_1,m_2,m_3,n_4=0}^{\infty}{\frac{m_2\left(m_2-1\right)}{\left(m_1+1\right)\left(m_1+2\right)}c_{m_1,m_2,m_3,m_4}x_1^{m_1+2}x_2^{m_2-2}x_3^{m_3}x_4^{m_4}}$$
Obviously, this formula only equals 0 if no values of $m_2$ other than 0 or 1 are permitted. The same logically applies for the other coordinates. As such, $f\left(\vec{x}\right)$ must have the form
$$f\left(\vec{x}\right)=f_{000}\left(x_1\right)+x_2\left(f_{100}\left(x_1\right)+x_3\left(f_{110}\left(x_1\right)+x_4f_{111}\left(x_1\right)\right)+x_4f_{101}\left(x_1\right)\right)+x_3\left(f_{010}\left(x_1\right)+x_4f_{011}\left(x_1\right)\right)+x_4f_{001}\left(x_1\right)$$
which cannot describe all metric tensor fields. Therefore, $\widetilde{g}=\nabla\vec{v}+\left(\nabla\vec{v}\right)^T$ doesn’t describe all metric tensor fields. These finding are relevant to all linear formulas. If they cannot produce a field for one unique component, they aren’t general. The search continues! (or ends here)
Best Answer
Both identities are correct. The point is that they refer to different bases. The former uses the decomposition, $$\vec{V} = V^r \partial_r + V^\theta \partial_\theta + V^z \partial_z\:,$$ the latter uses the decomposition, $$\vec{V} = V'^r \vec{e}_r + V'^\theta \vec{e}_\theta + V'^z \vec{e}_z\:,$$ where $\vec{e}_i = \frac{1}{\sqrt{g_{ii}}} \partial_{x^i}$ are unit vectors, so that $V'^i = \sqrt{g_{ii}} V^i$ (there is no sum over the repeated index $i$).