By symmetry:
$$I=4\int _0^1\int _0^1\int _0^x\int _0^z\sqrt{(z-w)^2+(x-y)^2} \, dw \, dy \, dx dz$$
Substituting $w \to z w, \quad y \to x y$:
$$I=4\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{z^2(1-w)^2+x^2(1-y)^2}~ x z \, dw \, dy \, dx dz$$
Substituting $w \to 1-w, \quad y \to 1-y$:
$$I=4\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{z^2 w^2+x^2 y^2}~ x z \, dw \, dy \, dx dz$$
Substituting $x^2 = u, \quad z^2 = v$:
$$I=\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{v w^2+u y^2}~ \, dw \, dy \, du dv$$
Substituting $v w^2=p, \quad u y^2=q$:
$$I=\int _0^1\int _0^1\int _0^{y^2}\int _0^{w^2} \sqrt{p+q}~ \, dp dq \frac{dy dw}{y^2 w^2}$$
$$I=\frac{2}{3} \int _0^1\int _0^1\int _0^{y^2}\left((q+w^2)^{3/2}-q^{3/2} \right) dq \frac{dy dw}{y^2 w^2}$$
$$I=\frac{4}{15} \int _0^1\int _0^1\left((y^2+w^2)^{5/2}-y^5-w^5 \right) \frac{dy dw}{y^2 w^2}$$
By symmetry:
$$I=\frac{8}{15} \int _0^1\int _0^w\left((y^2+w^2)^{5/2}-y^5-w^5 \right) \frac{dy dw}{y^2 w^2}$$
Substituting: $y = w s$:
$$I=\frac{8}{15} \int _0^1\int _0^1 w^2 \left((1+s^2)^{5/2}-s^5-1 \right) \frac{ds dw}{s^2 }$$
$$I=\frac{8}{45} \int _0^1 \left((1+s^2)^{5/2}-s^5-1 \right) \frac{ds }{s^2 }$$
The last single integral can be evaluated by different methods, and the result agrees with the one from the OP.
One good method here is integration by parts. I leave this as an exercise.
That $\cosh^2$ reminds me of one integral representation for the $\zeta$ function. Let's start from scratch: for any $s>0$,
$$ \sum_{n\geq 1}\frac{(-1)^{n+1}}{n^s} = \frac{1}{\Gamma(s)}\int_{0}^{+\infty}\frac{x^{s-1}}{e^x+1}\,dx $$
is a straightforward consequence of the (inverse) Laplace transform. If we apply integration by parts we get that
$$ \eta(s) = \frac{1}{\Gamma(s+1)}\int_{0}^{+\infty}\frac{x^s e^{x}}{(e^x+1)^2}\,dx = \frac{2^{s-1}}{\Gamma(s+1)}\int_{0}^{+\infty}\frac{x^s}{\cosh^2(x)}\,dx $$
holds for any $s>-1$. Similarly,
$$ \beta(s)=\sum_{n\geq 0}\frac{(-1)^n}{(2n+1)^s}=\frac{1}{2^{s+1}\Gamma(s)}\int_{0}^{+\infty}\frac{z^{s-1}}{\cosh(z/2)}\,dz $$
for any $s>0$ leads to
$$ \beta(s) = \frac{1}{2\Gamma(s+1)}\int_{0}^{+\infty}\frac{z^s\sinh(z) }{\cosh^2(z)}\,dz $$
for any $s>-1$. It follows that
$$ \frac{\sqrt{\pi}}{2}\beta\left(\tfrac{1}{2}\right) = \frac{1}{2}\int_{0}^{+\infty}\frac{\sqrt{z}\sinh(z)}{\cosh^2(z)}\,dz=\int_{0}^{+\infty}\frac{z^2\tanh(z^2)}{\cosh(z^{\color{red}{2}})}\,dz. $$
Long story short: there's a typo.
Best Answer
One can exploit a symmetry and split the integral in half along the line $y=x$ then use the following modified polar coordinates:
$$x = s^{\frac{2}{3}}\cos\theta \hspace{10 pt} y = s^{\frac{2}{3}}\sin\theta$$
$$\implies 2\cdot \frac{2}{3}\int_0^{\frac{\pi}{4}} \int_0^{\sec^{\frac{3}{2}}\theta} \frac{\log\left(s^{\frac{4}{3}}\right)}{\sqrt{\cos\theta+\sin\theta}}dsd\theta = \frac{16}{9}\int_0^{\frac{\pi}{4}} \frac{\sec^{\frac{3}{2}}\theta}{\sqrt{\cos\theta+\sin\theta}}\left(\frac{3}{2}\log(\sec\theta)-1\right)d\theta$$
then let $x=\tan\theta$
$$ \implies \frac{4}{3}\int_0^1 \frac{\log(1+x^2)}{\sqrt{1+x}}dx - \frac{16}{9}\int_0^1 \frac{1}{\sqrt{1+x}}dx$$
The integral on the right evaluates to $\frac{32}{9}(\sqrt{2}-1) \equiv \frac{32}{9}a$. The integral on the left becomes
$$ = \frac{8}{3}\sqrt{2}\log 2 - \frac{16}{3}\int_0^1 \frac{x\sqrt{1+x}}{1+x^2}dx$$
We have one and a half of the three terms the user posted, but this last integral is tricky and isn't yielding to a variety of methods. I will try to finish later, but in the meantime if anyone has any clever suggestions for this last integral I will be happy try it.