Why $r^{2} $ instead of $r^3$ in this change of variable?

What you have is a surface integral. In a surface integral you deal with a parametrization (in the case of your sphere) $$ \sigma(\theta,\phi)=(x_1+r\cos\theta\sin\phi,x_2+r\sin\theta\sin\phi,x_3+r\cos\phi),\ \ 0\leq\theta\leq2\pi,\ \ 0\leq\phi\leq\pi. $$ And your integral is calculated by $$\tag1 \int_{\partial B(x,r)}f(\sigma)\,d\sigma=\int_0^\pi\int_0^{2\pi}f(\sigma(\theta,\phi))\,\|\sigma_\theta\times\sigma_\phi\|\,d\theta\,d\phi. $$ You can rewrite $\sigma$ as $$ \sigma(\theta,\phi)=x+r\,\omega(\theta,\phi), $$ where $\omega$ is the parametrization of the unit ball centered at the origin. If you calculate $\omega_\theta\times\omega_\phi$, you'll see that $$ \sigma_\theta\times\sigma_\phi=r^2\,\omega_\theta\times\omega_\phi. $$ So $(1)$ becomes \begin{align} \int_0^\pi\int_0^{2\pi}f(\sigma(\theta,\phi))\,\|\sigma_\theta\times\sigma_\phi\|\,d\theta\,d\phi&=\int_0^\pi\int_0^{2\pi}f(x+r\omega(\theta,\phi))\,r^2\,\|\omega_\theta\times\omega_\phi\|\,d\theta\,d\phi\\ \ \\ &=\int_{\partial B(0,1)} f(x+r\omega)\,r^2\,d\omega. \end{align}

$\omega$ is confined to $\partial B(0,1)$. So the map you define $T$ is not from $\mathbb{R}^3$ to $\mathbb{R}^3$. Rather it is from $\partial B(0,1)$ to $\partial B(x,r)$. So there will be some $2\times 2$ Jacobian, depending on parametrization.

So far, this is just an attempt to explain why $r^3$ is not the determinant of the Jacobian.

So why does it actually work out to $r^2$? I will try to avoid some paramterization of these surfaces. Instead, recall that a determinant expresses the factor by which volume (or area in 2D, or hyper-volume...) is scaled. If you can picture your map $T$, it is uniformly scaling the sphere's surface area by a factor of $r^2$. Whatever parameterization you choose, calculating the determinant of the Jacobian needs to simplify down to $r^2$ with this map $T$.