A standard technique of mine when evaluating iterated integrals, such as double integrals, is to invoke a substitution that relates the integration variables of both integrals. To demonstrate and clarify what I mean, I shall present three examples for double integrals:
Example 1: Gaussian integral, $\Gamma\left(\frac{1}{2}\right)=\sqrt{\pi}$ $$I=\Gamma \left(\frac{1}{2}\right)^2 = \int_{0}^{\infty} \int_{0}^{\infty} \frac{1}{\sqrt{x t}}e^{-(t+x)}\,dx\,dt$$ We now enforce the substituion $t=sx \implies dt=x\,ds$: $$\implies I = \int_{0}^{\infty}\int_{0}^{\infty}\frac{1}{\sqrt{s}}e^{-x (s+1)}\,dx\,ds=\int_{0}^{\infty}\frac{1}{\sqrt{s}(s+1)} \,ds=\pi$$ $$\implies \Gamma\left(\frac{1}{2}\right)=\int_{-\infty}^{\infty}e^{-x^2}\,dx=\sqrt{\pi}$$ This is a much more elementary approach than the common Jacobian polar coordinates method of evaluating the Gaussian integral and is how I originally discovered the solution for myself a few years ago prior to learning about the Jacobian.
Example 2: In a previous answer of mine, I determined the following result: $$J=\frac{1}{2\pi}\int_{-\infty}^{\infty} \exp\left(\frac{x^2}{2}-\frac{y^2}{2}\right)K_0\left(\sqrt{(x-y)^2+\left(\frac{y^2}{2}-\frac{x^2}{2}\right)^2}\right)\,dy = \frac{e}{2} \sqrt{\pi}\operatorname{erfc}(1)$$ and is independent of $x$. Here $K_0$ is the modified Bessel function of the second kind of order $0$.
A crucial part of my method of evaluation was to invoke a particular substitution on the integral: $$J=\frac{1}{2\pi}\int_0^\infty\int_0^\infty \frac{\cosh\left(xv\left(\frac{v^2}{4t}+1\right)\right)}{t}\exp\left(-t-\frac{1+2t+x^2}{4t}v^2-\frac{v^4}{16t}\right)\,dv\,dt$$
We now enforce the substitution $t = s v^2 \implies dt = v^2 \,ds$: $$\begin{align}J&=\frac{1}{2\pi} \int_{0}^{\infty} \int_{0}^{\infty}\frac{\cosh\left(x v \left(\frac{1}{4s}+1\right)\right)}{s} \exp\left(-sv^2-\frac{1+2sv^2+x^2}{4s}-\frac{v^2}{16s}\right) \, dv \, ds\\&=\frac{1}{2\pi}\int_{0}^{\infty}\int_{0}^{\infty}\frac{\cosh\left(xv\left(\frac{1}{4s}+1\right)\right)}{s}\exp\left(-\frac{1+x^2}{4s}\right)\exp\left(-v^2 \frac{(1+4s)^2}{16s}\right)\,dv\,ds\end{align}$$ This allows the result to fall out readily, since: $$\int_{0}^{\infty}\cosh(a v)\exp\left(-v^2 b^2\right) \, dv=\frac{\sqrt{\pi}}{2b}\exp\left(\frac{a^2}{4b^2}\right)$$
$$\implies J = \frac{1}{2\pi} \int_{0}^{\infty} \frac{\exp \left(-\frac{1+x^2}{4s}\right)}{s} \cdot \frac{\sqrt{\pi}}{2\left(\frac{1+4s}{4\sqrt{s}}\right)}\exp \left(\frac{4s x^2 \left(\frac{1}{4s}+1\right)^2}{(1+4s)^2}\right)\, ds$$ This results in the integral being independent as we wanted since the $x^2$ terms cancel.
$$\implies J = \frac{1}{\sqrt{\pi}} \int_{0}^{\infty} \frac{\exp\left(-\frac{1}{4s}\right)}{\sqrt{s} (1+4s)} \, ds\stackrel{s\,\mapsto\frac{1}{s}}{=}\frac{1}{\sqrt{\pi}} \int_{0}^{\infty} \frac{\exp\left(-\frac{s}{4}\right)}{\sqrt{s} (4+s)} \, ds=\frac{e}{2} \sqrt{\pi} \, \operatorname{erfc} (1)$$
Example 3: Derivation of the relationship between the Beta function and Gamma function. $\Gamma(x)\Gamma(y)=\Gamma(x+y)B(x,y)$: $$\Gamma(x)\Gamma(y)=\int_{0}^{\infty}\int_{0}^{\infty} e^{-u-v}u^{x-1}v^{y-1}\,du\,dy$$ We now enforce the substitution $u=zt$ and $v=z(1-t)$: $$\begin{align}\Gamma(x)\Gamma(y)&=\int_{0}^{\infty}\int_{0}^{1}e^{-z}(zt)^{x-1}(z(1-t))^{y-1} z\,dt\,dz\\ &=\int_{0}^{\infty}e^{-z}z^{x+y-1}\,dz \cdot \int_{0}^{1}t^{x-1}(1-t)^{y-1}\\&=\Gamma(x+y)B(x,y)\end{align}$$ This is interestingly the method mentioned in the derivation of this property on the Wikipedia page of the Beta function. In this question about this derivation, the comments and answer invoke the Jacobian, suggesting that such substitutions can be made rigorous using the Jacobian, however, I would appreciate seeing exactly how that works as I only know about the Jacobian based on my own self-teaching.
So to summarise is it always valid to perform substitutions on iterated integrals, such as double or triple integrals in this way, whereby one or more integration variables are written in terms of an integration variable of a different integral in the iteration. I have never come across any issue in using this technique in the past, but have never found rigorous proof in any literature that states that such substitutions are valid; it surprises me that such a simple solution to the Gaussian integral is not often used.