I would like to compute the area of the sphere of radius $R$ with the following method:
I consider section of the sphere in an orthonormal frame with Cartesian coordinates $(x,y,z)$ delimited by the two horizontal planes of elevation $z$ and $z+dz$ so that the area of this section (which is a trapezoid) is $dS = dz/2\; (2\pi r(z)+2\pi r(z+dz))$ where $r(z)$ is the radius of the circle defined by the intersection of the Sphere and the plane of elevation $z$ given by $$ r(z)=\sqrt{R^2-z^2}. $$
So the integral to compute is $$ A = \pi\int_{-R}^R (r(z)+r(z+dz)) dz $$
I can write $r(z+dz)=r(z)+dz\,r'(z) + o(dz)$ and then
$$ A = \pi\int_{-R}^R (2r(z)+ dr +o(dz))dz $$
Now how to handle for instance $dr\,dz$ ? $dz\, dz$ ?
I don't master exterior derivative of differential forms, but how can we introduce them here ?
I can also do the following computation $$ A = \pi\int_{-R}^R (2r(z)+ r'(z)dz +o(dz))dz \approx \pi\int_{-R}^R (2r(z)+ r'(z)dz)dz $$ but here I have also to deal with $dz\, dz$. If it is equal to $dz\wedge dz=0$, I would like to see why geometrically.
The error in the analysis is in your approximation by the area of a trapezoid. The area of a frustum of a cone is given by the product of the average circumference with the slant height (as opposed to vertical height). [You can see the error in your approach if you try to do arclength ignoring the hypotenuse of the right triangle.]
You will need to do a surface integral to bring in double integration, which is not what you're trying to do.