I'm familiar with Riemann sums and using them to define single, double, and triple integrals in Cartesian coordinates. I also understand the idea that if I transform a coordinate system I will also be transforming the shape and size of the infinitesimal area elements that make up the infinitely fine grid over which we are doing our Riemann sums. So it makes intuitive sense to me that some transformation of the differentials would show up when we express integrals in one coordinate space in terms of integrals in a different coordinate space.
The part that I can't find an intuitive justification for is why multiplying by the absolute value of the Jacobian determinant specifically gives the area of the new infinitesimal element. That is why is it when taking double integrals:
$$dA = dxdy = \left|\begin{bmatrix}\frac{\partial x}{\partial u} & \frac{\partial y}{\partial u} \\ \frac{\partial x}{\partial v} & \frac{\partial y}{\partial v}\end{bmatrix}\right|dudv$$
Preferably I'm looking for an explanation that just justifies the formula that you get from expanding the determinant formula so that it can be equally well understood by people not familiar with linear-algebra/basis.
The determinant of a $2\times2$ matrix can be interpreted as the (signed) area of a parallelogram with sides defined by the columns or rows of the matrix. One way to see that this is so is to use some properties of determinants: by adding a suitable multiple of the first column to the second you can make them orthogonal without changing the value of the determinant, then apply a rotation to align the two columns with the coordinate axes, which also leave the determinant unchanged. The columns of the resulting diagonal matrix define an axis-aligned rectangle, and the values along the main diagonal are its side lengths. Now, if the matrix represents a linear transformation, it maps the unit square to a parallelogram, and you can view its determinant as the area-scaling factor of the transformation. More generally, the determinant of an $n\times n$ matrix is the signed volume of the polytope with edges defined by the columns/rows of the matrix.
Each term in the Riemann sum for a double integral consists of the integrand evaluated at some point on a rectangular mesh cell times the area of the cell. A coordinate transformation $\phi:(u,v)\mapsto(x,y)$ will in general map this cell to an area bounded by curved lines. The Jacobian matrix $J={\partial(x,y)\over\partial(u,v)}$ represents a linear map that best approximates the change in $\phi$ at a point. This map transforms the rectangular mesh cell with width $\Delta u$ and height $\Delta v$ into a parallelogram, and the area of this parallelogram is $\Delta u\Delta v \det(J(u,v))$, i.e., per the above discussion the area of the rectangular cell is scaled by the Jacobian determinant of $\phi$ evaluated at a vertex of the cell. If we don’t care about orientation, i.e., we’re computing an absolute double integral $\iint \sigma\,dx\,dy$ and not $\iint f\,dx\wedge dy$, then we don’t care about tracking orientation changes and use the absolute value of the Jacobian determinant.
Now for a rigorous proof, you’d of course need to show that the errors introduced by this approximation vanish quickly enough as the mesh size goes to zero and to generalize to more than two dimensions, but the basic idea remains the same.