How can one show the multiple integral $\int_{0}^{a} \int_{0}^{a} \frac{dx ~ dy}{(a^2+x^2+y^2)^{3/2}}= \frac{\pi}{6a}$ by hand?

Question

Mathematica gives the value of the multiple integral $$ \int_{0}^{a} \int_{0}^{a} \frac{dx ~ dy}{(a^2+x^2+y^2)^{3/2}}=\frac{\pi}{6a}. $$

See also this result from WolframAlpha:

The question is how to show this by hand?

This is posted in the spirit of sharing some experience and to see how others would solve it. I am posting a solution below (now moved to the duplicate target).