Thanks for the tip. I have solved the problem now for general $r_i$ and $k$. Yes, like you said, I had to calculate three integrals from the unique combinations of two radius lengths. This corresponds to a "corner" of the superellipsoid, so then you need to multiply by 8 to recover the full surface area. And the integration bounds do actually change, so part of my previous comment was incorrect.

That's brilliant. I got the same result for $r_1=r_2=r_3=1$. But I have trouble generalizing to arbitrary $r_i$. I restricted the domain to 1/4 of the cube's top face so that $u, v\geq0$ and we don't have to consider parity of $k$. Then you multiply by 24 to find the full area. I have that $\lambda= \left[ \left(\frac{u}{r_1} \right)^k+\left(\frac{v}{r_2}\right)^k+\left(\frac{z}{r_3}\right)^k)\right]^{-1/k}$ and $\gamma$ is unchanged. Integration region is the same but the lower bounds are 0. When I do this, the area does not match with the expected areas at $k=2$ and $k \to \infty$.

@BiswajitBanerjee That is a good point. I went back to Mathematica and checked the integrand, but I got the same result. So I'm confident the integrand is theoretically correct but it does seem strange that it appears it goes to 0 as $k$ increases. I also tried to integrate with Mathematica's NIntegrate. It matches with other methods at $k=2, 3$ but also fails at large $k$. But it gives a more informative error: "Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small."