Suppose, for argument's sake, that we have a function $f: \mathbb{R}^2 \to \mathbb{R}$ that we want to integrate over a circular surface. We can use polar coordinates and cut up the $r$ and $\theta$ intervals into $N$ and $M$ pieces respectively, and then use that:
$$ \begin{aligned} \hat S &= \sum_i^{N - 1}\sum_j^{M-1} f(r_i, \theta_j)\Delta S \\ &= \sum_i^{N - 1}\sum_j^{M-1} \frac{1}{2}f(r_i, \theta_j)(\theta_{j+1} - \theta_j)(r_{i+1}^2 - r_i^2) \\ &= \sum_i^{N-1}\sum_j^{M-1} \frac{1}{2}f(r_i, \theta_j)\Delta\theta((r_i + \Delta r)^2 - r_i^2) \\ &= \sum_i^{N-1}\sum_j^{M-1} \frac{1}{2}f(r_i, \theta_j)\Delta\theta(2r_i\Delta r + (\Delta r)^2) \\ &= \sum_i^{N-1}\sum_j^{M-1} f(r_i, \theta_i)r_i\Delta\theta\Delta r + \frac{1}{2}f(r_i, \theta_j)\Delta\theta(\Delta r)^2 \\ \end{aligned} $$
Breaking this up, we can take the limit w.r.t. $N$ and $M$ and get:
$$ \begin{aligned} S &= \lim_{N, M \to \infty} \Big[\sum_i^{N-1}\sum_j^{M-1} f(r_i, \theta_i)r_i\Delta\theta\Delta r\Big] + \lim_{N, M \to \infty}\Big[\sum_i^{N-1}\sum_j^{M-1}\frac{1}{2}f(r_i, \theta_j)\Delta\theta(\Delta r)^2\Big] \\ &= \int_0^R\int_0^{2\pi} f(r, \theta)rdrd\theta + \lim_{N, M \to \infty}\Big[\sum_i^{N-1}\sum_j^{M-1}\frac{1}{2}f(r_i, \theta_j)\Delta\theta(\Delta r)^2\Big] \end{aligned} $$
Now, by the standard definition of the surface integral, the result should be:
$$ S = \int_0^R\int_0^{2\pi} f(r, \theta)rdrd\theta $$
however that leaves the question: Why is it that
$$ \lim_{N, M \to \infty}\Big[\sum_i^{N-1}\sum_j^{M-1}\frac{1}{2}f(r_i, \theta_j)\Delta\theta(\Delta r)^2\Big] = 0 $$
And in general, how am I supposed to approach Riemann sums with non-linear, arbitrary $\Delta x$ terms, e.g. $(\Delta x)^2$, $\sqrt{\Delta x}$, etc.?
I think I have found a way to prove this, although it involves working with the definition of a Riemann integral directly. I would be extremely happy if someone showed a simpler way to do this (if it exists).
Since the above example is a double integral where $\theta$ basically doesn't matter, I'm going to prove this for the general case of a function $f(x)$ over an interval $(a, b)$. Here it goes.
We want to show that for all $\epsilon$ there exists a given partition $\Pi = \{x_1, x_2, ...\}$ such that for any refinement $\Pi_1$ of $\Pi$
$$ \Bigg\lvert \sum_i^{\Pi_1} f(x_i)(\Delta x_i)^2 \Bigg\rvert < \epsilon $$
A refinement is simply the original partition with more cuts. Therefore, we only need to show that splitting one sub-interval into two sub-intervals gives a smaller sum than the original. First, consider that if $\lvert f(x) \rvert < M$ for all $x\in (a, b)$, then
$$ \Bigg\lvert \sum_i^{\Pi} f(x_i)(\Delta x_i)^2 \Bigg\rvert \leq \sum_i^{\Pi} M(\Delta x_i)^2 $$
Now let us introduce another splitting point to $\Pi$, say $x_j$, between two other existing points $x_p < x_q$. The above sum changes to
$$ \sum_i^{\Pi_1} M(\Delta x_i)^2 = \sum_i^{\Pi}M(\Delta x_i)^2 - M(x_q - x_p)^2 + M(x_q - x_j)^2 + M(x_j - x_p)^2 $$
Setting $\alpha = x_q - x_p$, $\beta = x_q - x_j$ and $\gamma = x_j - x_p$, we know that:
$$ \alpha = \beta + \gamma \implies \alpha^2 > \beta^2 + \gamma^2 $$
Therefore, we can say that:
$$ \Bigg\lvert \sum_i^{\Pi_1}f(x_i)(\Delta x_i)^2 \Bigg\rvert \leq \sum_i^{\Pi_1} M(\Delta x_i)^2 < \sum_i^{\Pi} M(\Delta x_i)^2 $$
Now we show that any $\epsilon$ can be expressed as $\sum_i^{\Pi} M_1(\Delta x_i)^2$ given the right $\Pi$. Essentially, we have to solve the following two equations:
$$ \sum_i^\Pi \Delta x_i = b - a \\ \sum_i^\Pi (\Delta x_i)^2 = \frac{\epsilon}{M_1} $$
If we take $\Delta x_i = \frac{b-a}{N}$ for some $N$, then
$$ \sum_i^\Pi (\Delta x_i)^2 = N\frac{(b - a)^2}{N^2} = \frac{(b - a)^2}{N} = \frac{\epsilon}{M_1} $$
Rewriting this, we have
$$ \epsilon = \frac{M_1(b - a)^2}{N} $$
We need to show that we can choose a natural number $N$ and a real number $M_1$ such that $\lvert f(x)\rvert \leq M \leq M_1$ and the above is possible. If $\frac{\epsilon}{(b - a)^2} \geq M$, then taking $M_1 = \frac{\epsilon}{(b - a^2)}$ and $N = 1$ is sufficient. If $\frac{\epsilon}{(b - a)^2} < M$, then we take $N = \Big\lceil \frac{M(b-a)^2}{\epsilon} \Big\rceil$ and $M_1 = \frac{N\epsilon}{(b - a)^2}$, which gives us
$$ \epsilon = \frac{\frac{N\epsilon}{(b - a)^2}(b-a)^2}{N} = \epsilon $$
and
$$ M_1 = \frac{\lceil \frac{M(b-a)^2}{\epsilon}\rceil\epsilon}{(b-a)^2} \geq \frac{\frac{M(b-a)^2}{\epsilon}\epsilon}{(b-a)^2} = M \implies M_1 \geq M $$
Finally, we can say that for any $\epsilon$, we can construct a partition $\Pi$ such that:
$$ \sum_i^{\Pi}f(x_i)(\Delta x_i)^2 \leq \sum_i^{\Pi}M_1(\Delta x_i)^2 = \epsilon $$
Above we showed that whenever we refine a partition, the resulting sum is always less than the $M_1$-bound of the original partition, therefore for every refinement $\Pi_1$ of $\Pi$
$$ \Bigg\lvert \sum_i^{\Pi_1}f(x_i)(\Delta x_i)^2 \Bigg\rvert < \epsilon $$
Thus we conclude that the above sum converges to $0$, or in other words if the mesh of the partition is $P$
$$ \lim_{P\to 0} \sum_i^{\Pi} f(x_i)(\Delta x_i)^2 = 0 $$
Q.E.D.
On the question about having $\sqrt(\Delta x)$ etc., I believe that the above method would yield the same results.