$\newcommand{\d}{\,\mathrm{d}}$It can be shown that arclength, considered as a sum of increasingly fine partitions of the graph, approaches the integral formulation. However, I have only ever seen the surface area integral treated as a definition, or as a truth, and today I wish to prove it – my proof is quite unwieldy, so I'm not sure if it's rigorous or if it can be cleaned up. I also don't study differential geometry, so I'd like to know if my reasoning and terminology is consistent with the field. As soon as I saw the above formula I knew I had to attempt a proof, so here it is:
Definitions:
Let $M$ be a smooth surface if: $M\subset\Bbb R^3$ is bounded and of null Lebesgue measure and there exists a simply connected open $\Omega\subset\Bbb R^2$ together with a $C^2$ bijection $\varphi:\Omega\to M$.
Let $P=\{x_{i,j}\in M:i,j\in\{1,\cdots,N\}\}$ be a partition of a smooth surface $M$ if $\{x_{i,j}:j\in\{1,\cdots,N\}\}$ partitions a contour in $M$ for any fixed $1\le i\le N$ and visa versa. Define the area of $P$ to be the sum of areas of the triangles with vertices $\{x_{i,j},x_{i,j+1},x_{i+1,j}\}$ and $\{x_{i+1,j},x_{i+1,j+1},x_{i,j+1}\}$.
Note: I chose to define it via triangles as triangles are always coplanar, whereas when I attempted this before with quadrilaterals I ran into problems – decomposing the quadrilaterals into two triangles seems to work!
Define the surface area of a smooth surface $M$ to be the limit of the areas of the partitions of $M$ as $\max\{\max_{i,j}\{x_{i+1,j}-x_{i,j}\},\max_{i,j}\{x_{i,j+1}-x_{i,j}\}\}\to0$, if it exists.
If $S\subset\Bbb R^3$ can be decomposed as the disjoint union of finitely many $M_1,M_2,\cdots,M_n$ smooth surfaces with parametristations $(\Omega_1,\varphi_1),\cdots,(\Omega_n,\varphi_n)$, then it has a surface area equivalent to the sum of the surface areas of the $M_k$.
I am a touch worried that there may be pathological surfaces for which this doesn't work but I think this is all ok so far.
Let $Q$ be an arbitrary rectangular partition into points $(u_i,v_i)_{i=1}^N$ of a smooth surface $\Omega$, and put $\Delta Q=\max\{\max_i\{u_{i+1}-u_i\},\max_i\{v_{i+1}-v_i\}\}$. As $\varphi$ is a continuous bijection, $\varphi(Q)$ will be a valid partition of $M=\varphi(\Omega)$. Consider for any fixed $i$ the two triangles described by the points $(u_i,v_j),(u_{i+1},v_j),(u_i,v_{j+1}),(u_{i+1},v_{j+1})$. The corresponding area of the two triangles in the partition of $M$, $P=\varphi(Q)$, will be: $$\frac{1}{2}\|(\varphi(u_{i+1},v_j)-\varphi(u_i,v_j))\times(\varphi(u_i,v_{j+1})-\varphi(u_i,v_j))\|$$And: $$\frac{1}{2}\|(\varphi(u_{i+1},v_{j+1})-\varphi(u_i,v_{j+1}))\times(\varphi(u_{i+1},v_{j+1})-\varphi(u_{i+1},v_j))\|$$
EDIT: As per Matteo's comment, I've now realised that the MVT fails not only for functions with multidimensional domain, but also for functions $\Bbb R\to\Bbb R^n$. I think the bijectivity of $\varphi$ saves the proof, but it will require more work that I cannot straight away think of.
CONT:
As $\varphi$ is $C^2$, one may employ the mean value theorem twice by holding $u$ or $v$ fixed and letting the other vary, making the function defined $\Bbb R\to\Bbb R^3$, so defining $\Delta u_i=u_{i+1}-u_i$ and similarly for $\Delta v_j$ gives: $$\begin{align}\varphi(u_{i+1},v_{j+1})-\varphi(u_i,v_{j+1})&=\varphi(u_{i+1},v_j)-\varphi(u_i,v_j)+\varphi(u_{i+1},v_{j+1})\\&\hspace{10pt}-\varphi(u_{i+1},v_j)-(\varphi(u_i,v_{j+1})-\varphi(u_i,v_j))\\&=\varphi(u_{i+1},v_j)-\varphi(u_i,v_j)+\Delta v_j\cdot(\partial_v\varphi(u_{i+1},\theta_v)-\partial_v\varphi(u_i,\theta_v))\\&=\varphi(u_{i+1},v_j)-\varphi(u_i,v_j)+\Delta v_j\Delta u_i\cdot\partial^2_{uv}(\theta_u,\theta_v)\\\varphi(u_{i+1},v_{j+1})-\varphi(u_{i+1},v_j)&=\cdots\\&=\varphi(u_i,v_{j+1})-\varphi(u_i,v_j)+\Delta u_i\Delta v_j\cdot\partial^2_{vu}\varphi(\vartheta_u,\vartheta_v)\end{align}$$
For mean values $\theta_u,\vartheta_u\in(u_i,u_{i+1}),\,\theta_v,\vartheta_v\in(v_j,v_{j+1})$.
Then the second area becomes: $$\frac{1}{2}\|(\varphi(u_{i+1},v_j)-\varphi(u_i,v_j))\times(\varphi(u_i,v_{j+1})-\varphi(u_i,v_j))\|\\\hspace{80pt}\\\hspace{20pt}+\frac{1}{2}\Delta u_i\Delta v_j\left\{\|\partial^2_{uv}\varphi(\theta_u,\theta_v)\times(\varphi(u_i,v_{j+1})-\varphi(u_i,v_j))\|\\\hspace{40pt}+\|(\varphi(u_{i+1},v_j)-\varphi(u_i,v_j))\times\partial^2_{vu}\varphi(\vartheta_u,\vartheta_v)\|\\\hspace{30pt}+\Delta u_i\Delta v_j\cdot\|\partial^2_{uv}\varphi(\theta_u,\theta_v)\times\partial^2_{vu}\varphi(\vartheta_u,\vartheta_v)\|\right\}$$
Employ the mean value theorem once more with feeling to the expressions of the form $\varphi(u_i,v_{j+1})-\varphi(u_i,v_j)$, using mean values $\lambda_u,\lambda_v$ in the appropriate ranges, and the two areas become:
$$\frac{1}{2}\Delta u_i\Delta v_j\|\partial_u\varphi(\lambda_u,v_j)\times\partial_v\varphi(u_i,\lambda_v)\|$$
And:
$$\frac{1}{2}\Delta u_i\Delta v_j\|\partial_u\varphi(\lambda_u,v_j)\times\partial_v\varphi(u_i,\lambda_v)\|+\\\hspace{20pt}\frac{1}{2}\Delta u_i\Delta v_j\left\{\Delta v_j\cdot\|\partial^2_{uv}\varphi(\theta_u,\theta_v)\times\partial_v\varphi(u_i,\lambda_v)\|+\Delta u_i\cdot\|\partial_u\varphi(\lambda_u,v_j)\times\partial^2_{vu}\varphi(\vartheta_u,\vartheta_v)\|\\\hspace{30pt}+\Delta u_i\Delta v_j\cdot\|\partial^2_{uv}\varphi(\theta_u,\theta_v)\times\partial^2_{vu}\varphi(\vartheta_u,\vartheta_v)\|\right\}$$
The area of the partition $P=\varphi(Q)$ is then expressible as (by summing the two triangles per summand):
$$\sum_{i,j=1}^{N-1}\|\partial_u\varphi(\lambda_u,v_j)\times\partial_v(u_i,\lambda_v)\|\cdot\Delta u_i\Delta v_j\\+\frac{1}{2}\sum_{i,j=1}^{N-1}\left\{\Delta v_j\cdot\|\partial^2_{uv}\varphi(\theta_u,\theta_v)\times\partial_v\varphi(u_i,\lambda_v)\|+\Delta u_i\cdot\|\partial_u\varphi(\lambda_u,v_j)\times\partial^2_{vu}\varphi(\vartheta_u,\vartheta_v)\|\\\hspace{30pt}+\Delta u_i\Delta v_j\cdot\|\partial^2_{uv}\varphi(\theta_u,\theta_v)\times\partial^2_{vu}\varphi(\vartheta_u,\vartheta_v)\|\right\}\cdot\Delta u_i\Delta v_j$$
As $\Delta Q\to0$, the latter sum clearly goes to $0$ as it features quadratic $\Delta$ terms, and by considering the first sum as a Riemann sum over a double partition tagged with the values $\lambda$, in the passage to the limit (if it exists) we find that the surface area is:
$$A(M)=\iint_{\Omega}\left\|\frac{\partial\varphi}{\partial u}\times\frac{\partial\varphi}{\partial v}\right\|\d u\d v$$
And:
$$A(S)=\sum_{k=1}^n\iint_{\Omega_k}\left\|\frac{\partial\varphi_k}{\partial u}\times\frac{\partial\varphi_k}{\partial v}\right\|\d u\d v$$
Is this a correct proof? And is the preliminary definition a sensible one? Furthermore, I expect that the $C^2$ condition can be dropped but I don't see how. It would have certainly been easier to consider quadrilaterals rather than triangles, but my attempts at that failed on the realisation that four points need not be coplanar, and thus need not be a good approximation to the surface – finding the area and creating good partitions was tricky. Hopefully my triangle approach works!
Best Answer
Here is the famous Schwarz example of a triangular mesh inscribed in a finite cylinder and the total surface area goes to infinity as we move the planes of the hexagons closer and closer. I am not sure if this gives a counterexample to your procedure or not, but I'll let you ponder it.
P.S. The bijectivity of the mapping won't rescue the Mean Value Equality. Consider $\phi(t) = (\cos t,\sin t)$. Instead of the domain $[0,2\pi]$, take $[0,2\pi-\epsilon]$. Now $\|\phi(2\pi-\epsilon)-\phi(0)\|<\epsilon$, but $\|(2\pi-\epsilon)\phi'(\xi)\| = 2\pi-\epsilon$ for any $\xi\in (0,2\pi-\epsilon)$, so we cannot have $\phi(2\pi-\epsilon)-\phi(0) = (2\pi-\epsilon)\phi'(\xi)$ for any $\xi$.