Thanks to Prof. Pavel Grinfeld for pointing out the flaw in the derivation in the question (see comment here). The key is to distinguish between the normal to $\Gamma$ and the normal to $C$. (This type of distinction is also made between normals in the context of the divergence theorem applied to volumes on page 242 of Prof. Grinfeld's book).
Here is the normal, $\mathbf{n}$, to $\Gamma$.
\begin{equation}
\begin{aligned}
\mathbf{n} &= \frac{t^\alpha\mathbf{S}_\alpha}
{\sqrt{S_{\beta\gamma}t^\beta t^\gamma}}
\times \nu \\
&= \frac{t^\alpha}
{\sqrt{S_{\beta\gamma}t^\beta t^\gamma}}
\left(\mathbf{S}_\alpha \times \nu\right) \\
&= \frac{t^\alpha}
{\sqrt{S_{\beta\gamma}t^\beta t^\gamma}}
\epsilon_{\delta\alpha}\mathbf{S}^\delta \\
&= n_\delta\mathbf{S}^\delta
\end{aligned}
\end{equation}
Thus,
\begin{equation}
n_\alpha = \frac{\epsilon_{\alpha\delta}t^\delta}
{\sqrt{S_{\beta\gamma}t^\beta t^\gamma}}
\end{equation}
By analogy, the coefficients of the normal, $\bar{\mathbf{n}}$, to $C$, are given by
\begin{equation}
\bar{n}_\alpha = \frac{\bar{\epsilon}_{\alpha\delta}t^\delta}
{\sqrt{\bar{S}_{\beta\gamma}t^\beta t^\gamma}}
\end{equation}
Now, $\bar{\epsilon}_{\alpha\delta} = \frac{\epsilon_{\alpha\delta}}{\sqrt{S}}$
and $\bar{S}_{\beta\gamma} = \delta_{\beta\gamma}$, so that
\begin{equation}
\bar{n}_\alpha = \frac{1}{\sqrt{S}}
\frac{\sqrt{S_{\beta\gamma}t^\beta t^\gamma}}
{\sqrt{\delta_{\mu\upsilon}t^\mu t^\upsilon}}
n_\alpha
\end{equation}
Now, the derivation in the question looks like this:
\begin{equation}
\begin{aligned}
\int\limits_\Omega \nabla_\alpha v^\alpha \text{d}\Omega
&= \int\limits_A \frac{1}{\sqrt{S}}\frac{\partial}{\partial S^\alpha}
\left(v^\alpha\sqrt{S}\right) \sqrt{S}\;\text{d}A
\quad \text{by the Voss-Weyl formula} \\
&= \int\limits_A \frac{\partial}{\partial S^\alpha}
\left(v^\alpha\sqrt{S}\right) \text{d}A \\
&= \oint\limits_C v^\alpha \bar{n}_\alpha \sqrt{S}\; \text{d}C
\quad \text{by the divergence theorem in $\mathbb{R}^2$} \\
&= \oint\limits_C v^\alpha n_\alpha
\frac{\sqrt{S_{\beta\gamma}t^\beta t^\gamma}}
{\sqrt{\delta_{\mu\upsilon}t^\mu t^\upsilon}}\; \text{d}C \\
&= \oint\limits_\Gamma v^\alpha n_\alpha \text{d}\Gamma
\end{aligned}
\end{equation}
and all is well!
That, from my perspective, resolves the conflict alluded to in my question.
However, I also want to share another calculation, which does not use the tensor notation formally, that guided me in the above derivation (I seem to have to switch between the two ways of thinking, each guiding the other). The following is a demonstration of the divergence theorem on the surface $\Omega$, knowing the divergence theorem in the coordinate space.
Let $\hat{\Xi}$ be the coordinate map, $\tilde{S}$ be a parametrization of $C$, and $\tilde{\Xi}$, the corresponding parametrization of $\Gamma$ as shown in the figure.
\begin{equation}
\begin{aligned}
& \oint_\limits\Gamma \mathbf{v}\cdot\mathbf{n}\text{d}\Gamma \\
= & \int_{t_0}^{t_1} \mathbf{v}\cdot
\left(\frac{\tilde{\Xi}'(t)}
{\left\Vert\tilde{\Xi}'(t)\right\Vert}
\times\nu \right)
\left\Vert\tilde{\Xi}'(t)\right\Vert \text{d}t \\
= & \int_{t_0}^{t_1} \mathbf{v}\cdot
\left(\tilde{\Xi}'(t) \times\nu \right) \text{d}t \\
= & \int_{t_0}^{t_1} \mathbf{v}\cdot
\left(
\left[\begin{array}{c|c}
\frac{\partial\hat{\Xi}}{\partial S^1} &
\frac{\partial\hat{\Xi}}{\partial S^2}
\end{array}\right] \tilde{S}'(t)
\times\nu
\right) \text{d}t \\
= & \int_{t_0}^{t_1} \mathbf{v}\cdot
\left(
\left[\begin{array}{c|c}
\frac{\partial\hat{\Xi}}{\partial S^1} \times \nu &
\frac{\partial\hat{\Xi}}{\partial S^2} \times \nu
\end{array}\right] \tilde{S}'(t)
\right) \text{d}t \\
= & \int_{t_0}^{t_1} \left(
\left[\begin{array}{c|c}
\frac{\partial\hat{\Xi}}{\partial S^1} &
\frac{\partial\hat{\Xi}}{\partial S^2}
\end{array}\right]
\begin{pmatrix} v^1 \\ v^2 \end{pmatrix}
\right)
\cdot
\left(
\left[\begin{array}{c|c}
\frac{\partial\hat{\Xi}}{\partial S^1} \times \nu &
\frac{\partial\hat{\Xi}}{\partial S^2} \times \nu
\end{array}\right] \tilde{S}'(t)
\right) \text{d}t \\
= & \int_{t_0}^{t_1} \begin{pmatrix} v^1 \\ v^2 \end{pmatrix}^\top
\left[\begin{array}{c|c}
\frac{\partial\hat{\Xi}}{\partial S^1} &
\frac{\partial\hat{\Xi}}{\partial S^2}
\end{array}\right]^\top
\left[\begin{array}{c|c}
\frac{\partial\hat{\Xi}}{\partial S^1} \times \nu &
\frac{\partial\hat{\Xi}}{\partial S^2} \times \nu
\end{array}\right] \tilde{S}'(t)
\text{d}t \\
& \quad \text{the integrand here can be recognized as}\
(v^1\mathbf{S}_1 + v^2\mathbf{S}_2) \cdot
\left[-\mathbf{S}^2 \; \mathbf{S}^1\right]
\tilde{S}'(t) S \\
= & \int_{t_0}^{t_1} \begin{pmatrix} v^1 \\ v^2 \end{pmatrix}^\top
\left[\begin{array}{cc}
0 & S \\
-S & 0
\end{array}\right] \tilde{S}'(t)
\text{d}t \\
= & \int_{t_0}^{t_1} \begin{pmatrix} v^1 \\ v^2 \end{pmatrix}^\top
\begin{pmatrix}
\tilde{S}^2{}'(t) \\ -\tilde{S}^1{}'(t)
\end{pmatrix}
S \text{d}t \\
= & \int_{t_0}^{t_1} \begin{pmatrix} v^1 S \\ v^2 S \end{pmatrix}^\top
\frac{
\begin{pmatrix}
\tilde{S}^2{}'(t) \\ -\tilde{S}^1{}'(t)
\end{pmatrix}
}{\left\Vert \tilde{S}'(t) \right\Vert}
\left\Vert \tilde{S}'(t) \right\Vert \text{d}t \\
= & \oint_\limits{C} \begin{pmatrix} v^1 S \\ v^2 S \end{pmatrix}^\top
\bar{\mathbf{n}} \text{d}C \\
= & \int_\limits{A}\left(
\frac{\partial}{\partial S^1}\left(v^1 S\right) +
\frac{\partial}{\partial S^2}\left(v^2 S\right)
\right)\text{d}A
\quad \text{by the divergence theorem in}\ \mathbb{R}^2 \\
= & \int_\limits{A}\frac{1}{S}
\left(
\frac{\partial}{\partial S^1}\left(v^1 S\right) +
\frac{\partial}{\partial S^2}\left(v^2 S\right)
\right) S \text{d}A \\
= & \int_\limits{\Omega}\frac{1}{S}
\left(
\frac{\partial}{\partial S^1}\left(v^1 S\right) +
\frac{\partial}{\partial S^2}\left(v^2 S\right)
\right) \text{d}\Omega
\end{aligned}
\end{equation}
which is the desired result.
Best Answer
First, you have to understand that not all bases are coordinate bases. A coordinate basis is one where the basis vectors are of the form $e_\mu = \partial_\mu = \frac{\partial}{\partial x^\mu}$ for some coordinate system $\{x^\mu\}_{\mu = 1}^d$. Equivalently, a coordinate basis is one in which the commutator $[e_\mu, e_\nu]$ between any pair of basis vectors vanishes. In differential geometry, ON bases (in which the metric is $g_{\mu\nu} = \delta_{\mu\nu}$) have a tendency not to be coordinate bases!
The regular polar coordinate system $(r,\theta)$ gives you a basis $\{\partial_r, \partial_\theta\}$ and a metric $g = \begin{pmatrix} 1 & 0 \\ 0 & r^2\end{pmatrix}$ with a determinant $\sqrt{\det(g)} = r$. In this basis, you would express a vector as $V = V^\mu \partial_\mu = V^r \partial_r + V^\theta \partial_\theta$.
Since this is a coordinate basis, we can use the given formula to evaluate the divergence: $$ \nabla_aV^a = \nabla_\mu V^\mu = \frac{1}{\sqrt{\det(g)}} \partial_\mu \left(\sqrt{\det(g)} V^\mu\right) = \frac{1}{r} \partial_\mu (rV^\mu) = \underline{\partial_r V^r + \frac{1}{r} V^r + \partial_\theta V^\theta}. $$
If you switch your basis to $\{e_1 = \partial_r, e_2 = \frac{1}{r} \partial_\theta\}$, the metric becomes $g_{\mu\nu} = \delta_{\mu\nu}$ and a vector is expressed as $V = V^\mu e_\mu = V^1 e_1 + V^2 e_2$, from which you can see that the new vector components are $V^1 = V^r$ and $V^2 = r V^\theta$. The new basis is not a coordinate basis; you can easily check that $[e_1, e_2] = -\frac{1}{r^2} \partial_\theta \neq 0$.
The expression $(\nabla_\mu V)^\nu = \partial_\mu V^\nu + \Gamma^\nu_{\rho\mu} V^\rho$ only holds in a coordinate basis. The symbol $\partial_\mu$ has no meaning in a non-coordinate basis. Instead, the general formula is $$ (\nabla_\mu V)^\nu = e_\mu(V^\nu) + \Gamma^\nu_{\rho\mu} V^\rho. $$ Here, $e_\mu(V^\nu)$ means that the differential operator $e_\mu$ acts on the function $V^\nu$. We cannot use the given formula for the divergence in this basis, because it only makes sense in a coordinate basis. Furthermore, even though the metric is $\delta_{\mu\nu}$ in this basis, the connection components $\Gamma^\mu_{\nu\rho}$ do not vanish since their expression as Christoffel symbols $\Gamma^\mu_{\nu\rho} = \frac{1}{2} g^{\mu\sigma} (g_{\nu\sigma,\rho} + g_{\sigma\rho,\nu} - g_{\nu\rho,\sigma})$ only makes sense in a coordinate basis.
To find the connection components in a general basis, we must use the general expression for the Levi-Civita connection, $$ Z \cdot \nabla_X Y = \frac{1}{2}\Big(X(Y \cdot Z) + Y(X \cdot Z) - Z(X \cdot Y) - X \cdot [Y,Z] - Y \cdot [X,Z] + Z \cdot [X,Y]\Big). $$ The connection components are defined as $\Gamma^\mu_{\nu\rho} = (\nabla_\rho e_\nu)^\mu = g^{\mu\sigma} (\nabla_\rho e_\nu)_\mu = g^{\mu\sigma} e_\mu \cdot \nabla_{e_\rho} e_\nu$, so we find (with $X = e_\rho$, $Y = e_\nu$, $Z = e_\sigma$) $$ \Gamma^\mu_{\nu\rho} = \frac{1}{2} g^{\mu\sigma} \Big( e_\rho (\overbrace{e_\nu \cdot e_\sigma}^{g_{\nu\sigma}}) + e_\nu (\overbrace{e_\rho \cdot e_\sigma}^{g_{\rho\sigma}}) - e_\sigma (\overbrace{e_\nu \cdot e_\rho}^{g_{\nu\rho}}) - e_\rho \cdot [e_\nu, e_\sigma] - e_\nu \cdot [e_\rho, e_\sigma] + e_\sigma \cdot [e_\nu, e_\rho] \Big). $$
Let us use this formula to compute $\nabla_\mu V^\mu$ in the basis $\{e_1, e_2\}$. Since the metric components in our case are constant, the first three terms vanish (these are the terms that give the Christoffel symbols in a coordinate basis). But the other three terms do not vanish. Inserting $g^{\mu\sigma} = \delta^{\mu\sigma}$, we get $$ \Gamma^\mu_{\nu\rho} = -\frac{1}{2} \Big( e_\rho \cdot [e_\nu, e_\mu] + e_\nu \cdot [e_\rho, e_\mu] - e_\mu \cdot [e_\nu, e_\rho] \Big). $$ As we found before, we have the commutator $[e_1, e_2] = -\frac{1}{r^2} \partial_\theta = -\frac{1}{r} e_2$. To compute $\nabla_\mu V^\mu = e_\mu(V^\mu) + \Gamma^\mu_{\nu\mu} V^\nu$, we need $\Gamma^\mu_{1\mu}$ and $\Gamma^\mu_{2\mu}$. Clearly $\Gamma^\mu_{\mu\mu}$ (no sum) vanishes because all commutators vanish, so $\Gamma^\mu_{1\mu} = \Gamma^2_{12}$ and $\Gamma^\mu_{2\mu} = \Gamma^1_{21}$. Doing the computation: $$ \begin{align} \Gamma^2_{12} &= -\frac{1}{2} \Big(e_2 \cdot [e_1, e_2] + e_1 \cdot [e_2, e_2] - e_2 \cdot [e_2, e_1]\Big) = -e_2 \cdot [e_1, e_2] = \frac{1}{r} e_2 \cdot e_2 = \frac{1}{r}, \\ \Gamma^1_{21} &= -\frac{1}{2} \Big(e_1 \cdot [e_2, e_1] + e_2 \cdot [e_1, e_1] - e_1 \cdot [e_1, e_2]\Big) = -e_1 \cdot [e_2, e_1] = -\frac{1}{r} e_1 \cdot e_2 = 0. \end{align} $$ We are then able to calculate the divergence: $$ \nabla_\mu V^\mu = e_1(V^1) + \overbrace{\Gamma^1_{21}}^0 V^2 + e_2(V^2) + \overbrace{\Gamma^2_{12}}^{1/r} V^1 $$ and inserting the expressions for $e_1$, $e_2$, $V^1$ and $V^2$: $$ \begin{align} \nabla_\mu V^\mu &= \partial_r(V^1) + \frac{1}{r} V^1 + \frac{1}{r} \partial_\theta(V^2) \\ &= \partial_r (V^r) + \frac{1}{r} V^r + \frac{1}{r} \partial_\theta (r V^\theta) \\ &= \underline{\partial_r V^r + \frac{1}{r} V^r + \partial_\theta V^\theta}, \end{align} $$ which is the same as before.
It is, in fact, always possible to find an orthonormal coordinate basis, at least locally. This is what is called a local inertial frame in GR: Starting from an ON basis $\{e_\mu\}_{\mu = 1}^d$ at a point $p$, you can define coordinates (normal coordinates) around $p$ via the geodesics through $p$. Starting from the basis $\{\partial_r, \partial_\theta\}$, we would get the regular Cartesian coordinate system, but centred on our starting point $p$, given by $(r_0, \theta_0)$: $$ \begin{align} \tilde{x} &= (r - r_0) \cos(\theta - \theta_0) \\ \tilde{y} &= (r - r_0) \sin(\theta - \theta_0). \end{align} $$ In this coordinate system, you will indeed find a diagonal metric as well as vanishing Christoffel symbols (generally only true locally, but here globally), which means that $\nabla_\mu V^\mu = \partial_\mu V^ \mu$.