Just fleshing out what @Ted Shifrin wrote in the comments. The case where the boundary is $C^1$ is straightforward, once you have the technical background (existence of partitions of unity) and you know the definition of what it means to have a $C^1$ boundary. However, if you want to avoid the language of differential forms, the argument appears to be more complicated because we have to manually keep track of various factors of determinants and so on (something which would have been automatic if we used differential forms). One could probably streamline to proof by more efficient use of charts on the boundary, but I didn't want to invoke too much analysis results regarding straightening of the boundary and canonical forms for vector fields on boundaries.
Note that this is all proven in Loomis and Sternberg's Advanced Calculus (for the divergence theorem they do things just an $\epsilon$ more generally, using densities). Pretty much the same proof is found in any differential geometry textbook for Stokes theorem; here I'm just rewording it to fit the divergence theorem.
Here's what we shall prove:
Let $(M,g)$ be an $n$-dimensional Riemannian manifold (say $n\geq 2$), of class $C^2$ and $U\subset M$ an open set with $C^1$ boundary $\partial U$, having continuous unit outward normal vector field $\nu$. If $X$ is a $C^1$ vector field on $M$ such that $\text{supp}(X)\cap \overline{U}$ is compact then we have
\begin{align}
\int_U\text{div}(X)\,dV&=\int_{\partial U} g(X,\nu)\, dS
\end{align}
where $dV=dV_g$ is the Riemannian measure on $M$ induced by the metric tensor $g$ and $dS=dV_{\iota^*g}$ is the Riemannian measure on $\partial U$ induced by the pull-back Riemannian metric $\iota^*g$ ($\iota:\partial U\to M$ being the $C^1$ canonical inclusion).
Step 1: Reduction to the Local Case.
Explicitly, since $\text{supp}(X)\cap \overline{U}$ is compact, we can find finitely many open sets $W_1,\dots, W_N$ in $M$, each having compact closure such that each $W_i$ is the domain of a $C^1$ chart map $\phi_i:W_i\to \phi_i[W_i]\subset\Bbb{R}^n$, and such that the $W_i$ fall into two categories:
- $W_i\subset U$ (meaning it is a chart for a point of $U$)
- $W_i$ contains a boundary point $p_i$, in which case the definition of $\partial U$ having a $C^1$ boundary tells us that we can choose $\phi_i[W_i]=(-1,1)^n$ and $\phi_i[W_i\cap U]= (-1,1)^{n-1}\times (0,1)$ (i.e the part containing $U$ lies in the upper half-space $\{x^n>0\}$).
Let $K=\bigcup_{i=1}^N\overline{W_i}$ and $W_0$ its complement in $M$; this is a compact set, and now $\{W_0,W_1,\dots, W_N\}$ is an open cover for $M$. We can thus find a ($C^2$ if you wish) partition of unity $\{h_0,h_1,\dots, h_N\}$ of $M$ subordinate to the above defined open cover. Since $\overline{U}\cap W_0=\emptyset$, it follows that $\sum_{i=1}^Nh_i=1$ on $\overline{U}$, so
\begin{align}
\int_U\text{div}(X)\,dV&=\int_U\text{div}\left(\sum_{i=1}^Nh_iX\right)\,dV=
\sum_{i=1}^N\int_U\text{div}(h_iX)\,dV
\end{align}
and similarly
\begin{align}
\int_{\partial U}g(X,\nu)\,dS&=\sum_{i=1}^N\int_{\partial U}g(h_iX,\nu)\,dS
\end{align}
Hence, it suffices to prove the divergence theorem for each of the $N$ $C^1$ vector fields $h_iX$. Notice that these vector fields fall into two categories, which we now discuss.
Step 2: Dealing with first type of vector fields.
The first case is having a $C^1$ vector field $Y$ on $M$, such that $\text{supp}(Y)$ is compact and lies inside the domain of a $C^1$ chart $(W,\phi=(x^1,\dots, x^n))$ of $M$, with $W\subset U$. In this case, the RHS of the divergence theorem is clearly $0$. The LHS is also easily seen to be $0$, since
\begin{align}
\int_U\text{div}(Y)\,dV=\int_W\text{div}(Y)\,dV&=\int_{\phi[W]}\sum_{i=1}^n\frac{1}{\sqrt{|g|}}\frac{\partial \left(\sqrt{|g|}Y^i\right)}{\partial x^i}\,\cdot \sqrt{|g|}\,dx^1\cdots dx^n\\
&= \int_{\phi[W]}\sum_{i=1}^n\frac{\partial \left(\sqrt{|g|}Y^i\right)}{\partial x^i}\,dx^1\cdots dx^n
\end{align}
here, I used the Voss-Weyl formula for the divergence of a vector field. Now each integral vanishes because $\sqrt{|g|}Y^i$ has compact support in $\phi[W]$; so for the $i^{th}$ term, use FUbini to integrate over $dx^i$ first, and then use the FTC (integral of derivative is difference at endpoints), and finally using compact support, it follows that these "endpoint" values are just $0$. Thus, for such a vector field, the divergence theorem holds because both sides are $0$.
Step 3: Dealing with the second type of vector fields.
Here, we're dealing with $C^1$ vector fields $Z$ on $M$ having compact support, which lies inside the domain of a $C^1$ chart $(W,\phi=(x^1,\dots, x^n))$, such that $\phi[W]=(-1,1)^{n}$, and $\phi[W\cap U]=(-1,1)^{n-1}\times (0,1)$. Very similarly to step 2, we have
\begin{align}
\int_U\text{div}(Z)\,dV&=\int_{W\cap U}\text{div}(Z)\,dV=\sum_{i=1}^n\int_{(-1,1)^{n-1}\times (0,1)}\frac{\partial \left(\sqrt{|g|}Z^i\right)}{\partial x^i}\,dx^1\cdots \,dx^n
\end{align}
For $i=1,\dots, n-1$ the integrals vanish for the exact same reason as in step 2: use Fubini to integrate over $dx^i$ from $(-1,1)$, and use the FTC, and compact support to see things vanish. For the $i=n$ term, we integrate over $dx^n$, but the crucial distinction is that we integrate over $x^n\in (0,1)$, not $(-1,1)$ (which is what happened in step 2, and hence this term also didn't contribute anything previously). Again, using Fubini, and FTC, the function vanishes as $x^n\to 1^-$ due to compact support, so we're only left with the contribution at $x^n=0$,
\begin{align}
\int_U\text{div}(Z)\,dV&=\int_{(-1,1)^{n-1}}-\left(\sqrt{|g|}Z^n\right)(x^1,\dots, x^{n-1},0)\,dx^1\dots\,dx^{n-1}.
\end{align}
(the minus sign is from the FTC). On the other hand, we have
\begin{align}
\int_{\partial U}g(Z,\nu)\,dS&=\int_{W\cap \partial U}g(Z,\nu)\,dS\\
&=\int_{(-1,1)^{n-1}}g_{ij}Z^i\nu^j\sqrt{\left|\det(g_{ab})_{1\leq a,b\leq n-1}\right|}\,dx^1\dots, dx^{n-1}\tag{$*$}
\end{align}
The whole integrand should be evaluated at $(x^1,\dots, x^{n-1},0)$, since that's what the boundary chart tells us to do. Now, note that $\phi$ maps $U\cap W$ into the upper-half plane described by $\{x^n>0\}$, and it maps $\partial U\cap W$ to $\{x^n=0\}$. Hopefully you recall that we find normal vectors to level sets by taking gradients of the functions defining the level set:
\begin{align}
\text{grad}_g(x^n)&= g^{ij}\frac{\partial x^n}{\partial x^i}\frac{\partial}{\partial x^j}=g^{nj}\frac{\partial}{\partial x^j}
\end{align}
and this vector field has squared norm equal to
\begin{align}
g\left(g^{nj}\frac{\partial}{\partial x^j}, g^{nk}\frac{\partial}{\partial x^k}\right)=g^{nj}g^{nk}g_{jk}=g^{nj}\delta^n_j=g^{nn}
\end{align}
Here, $[g^{ij}]$ is the inverse matrix of $[g_{ij}]$, so $g^{nn}$ refers to the $(n,n)$ entry of the inverse matrix of $[g_{ij}]$. Therefore, (evaluate at points of the boundary, i.e $(x^1,\dots, x^{n-1},0)$)
\begin{align}
\nu&=\frac{-\text{grad}_g(x^n)}{\|\text{grad}_g(x^n)\|}=-\frac{g^{nj}}{\sqrt{g^{nn}}}\frac{\partial}{\partial x^j}\tag{$**$}
\end{align}
Looking back at the integrand of $(*)$, we can plug in $(**)$ to get
\begin{align}
g_{ij}Z^j\nu^j\sqrt{\det (g_{ab})_{1\leq a,b\leq n-1}}&=g_{ij}Z^i\left(\frac{-g^{nj}}{\sqrt{g^{nn}}}\right)\sqrt{\det (g_{ab})_{1\leq a,b\leq n-1}}\\
&=-(g_{ij}g^{nj}Z^i)\sqrt{\frac{\det (g_{ab})_{1\leq a,b\leq n-1}}{g^{nn}}}\\
&=-(\delta_i^nZ^i)\sqrt{\det (g_{1j})_{1\leq i,j\leq n}}\\
&\equiv -Z^n\sqrt{|g|},
\end{align}
where we have used Cramer's formula for calculating inverse matrices: the $ij$ entry of the matri $A^{-1}$ is $\frac{1}{\det A}$ times the determinant of the matrix obtained by removing row $j$, column $i$ from $A$. This completes the proof of the equality in this case, and hence completes the entire proof of the divergence theorem.
Some remarks:
The benefit of using partitions of unity is that rather than chopping up the domain $U$ into a bunch of more manageable pieces (such as convex, or bounded by graphs of smooth functions, or any other "basic region" as you've seen in the other post), we are instead chopping up our vector field (or differential forms in Stokes' theorem) from $X$ to $\sum_{i=1}^Nh_iX$, so that each $h_iX$ is easily managed (i.e has support lying inside of a single nice chart).
Also, as I mentioned in the beginning, this proof of the divergence theorem appears long because I tried to avoid results about how nice things can look at the boundary, and also because I avoided the differential forms language. As a result, the difficulty had to be pushed elsewhere, and in this case, it is in explicitly calculating the surface integral in local coordinates (so the difficulty was pushed from analysis to linear algebra). The only real facts I used are Voss-Weyl formula for divergence and the local coordinate formula for the gradient of a function (and also that manifolds described locally as level sets have tangent space equal to the kernel of the function defining the level set, hence the gradient is orthogonal). In those links the manifolds are stated to be $C^{\infty}$, but we only need $C^1$ (after all we're only differentiating things once).
Along these same lines one can also prove Stokes theorem for a $C^2$ oriented $n$-dimensional manifold $M$, with $U$ an open set having $C^1$ boundary (with the induced boundary orientation), and a $C^1$ differential $(n-1)$-form $\omega$ on $M$ such that $\text{supp}(\omega)\cap \overline{U}$ is compact, in which case we still have $\int_Ud\omega=\int_{\partial U}\iota^*\omega$, where $\iota:\partial U\to M$ is the canonical injection, or simply written $\int_Ud\omega=\int_{\partial U}\omega$.
On an oriented $C^2$ Riemannian manifold $(M,g)$, the divergence theorem is indeed a special case of Stokes' theorem. If we let this time $dV,dS$ be the volume forms (not measures as before) on $M$ and $\partial U$, then For a $C^1$ vector field $X$, one can show that $\text{div}(X)$ satisfies
\begin{align}
\text{div}(X)\,dV&=L_X(dV)=d(\iota_XdV)+\iota_X(d\,dV)=d(\iota_XdV)
\end{align}
Apologies for the notation, but as always, $dV$ is not an exact form, but it does have top-degree, hence its exterior derivative vanishes (we also used Cartan's magic formula here). Hence, by Stokes theorem, we have
\begin{align}
\int_U\text{div}(X)\,dV&=\int_Ud(\iota_X\,dV)=\int_{\partial U}\iota_X(dV)=\int_{\partial U}\langle X,\nu\rangle_g\,dS,
\end{align}
where the final equality uses the stuff I explained in this answer
