How to find this four-dimensional integral as an iterated integral with Mathematica?

I mean the following integral $$\int\limits _{[0,1]^4}\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}\,dx_1dx_2dy_1dy_2 .$$ Its value is the average length of a random interval in the unit square (see, for example, MathWorld). Mathematica finds a multi-dimensional integral as an iterated integral. There are $4!=24$ variants to do it. Let us consider the first four items:

p = Permutations[{{x1, 0, 1}, {x2, 0, 1}, {y1, 0, 1}, {y2, 0, 1}}][[1 ;; 4]]; ParallelTable[Integrate[Sqrt[(x2 - x1)^2 + (y2 - y1)^2],i[[1]], i[[2]], i[[3]], i[[4]]], {i, p}] 

{1/240 (32 + 16 Sqrt[2] + (320 + 15 I) \[Pi] + 230 ArcSinh[1] + 30 ArcTanh[Sqrt[2]] + 30 Log[99 - 70 Sqrt[2]]), 1/240 (32 + 16 Sqrt[2] + (320 + 15 I) \[Pi] + 230 ArcSinh[1] + 30 ArcTanh[Sqrt[2]] + 30 Log[99 - 70 Sqrt[2]]), 1/240 (32 + 16 Sqrt[2] + (120 + 15 I) \[Pi] - 60 ArcCoth[Sqrt[2]] + 260 ArcSinh[1] + 30 ArcTanh[Sqrt[2]] + 25 Log[99 - 70 Sqrt[2]]), 1/240 (32 + 16 Sqrt[2] + (120 + 15 I) \[Pi] - 60 ArcCoth[Sqrt[2]] + 230 ArcSinh[1] + 30 ArcTanh[Sqrt[2]] + 20 Log[99 - 70 Sqrt[2]])}

N[%] 

{4.7102 + 0. I, 4.7102 + 0. I, 2.0922 + 0. I, 2.0922 + 0. I}

The calculation of each integral in 14.1 on Windows takes approximately 15 minutes for me. The rest of the iterated integrals produce similar results, some of those numerically equal 4.7102 + 0. I and some of those numerically equal 2.0922 + 0. I (the executed *.nb file on demand). All the results of the iterated integrations are not correct since the distance between two points in the unit square is less than or equal to $\sqrt 2$ and

NIntegrate[Sqrt[(x2-x1)^2+(y2-y1)^2],{x1,0,1},{x2,0,1},{y1,0,1},{y2,0,1},AccuracyGoal->8,PrecisionGoal->8,Method->"LocalAdaptive"] 

0.521404

and

Mean[TransformedDistribution[Sqrt[(x2 - x1)^2 + (y2 - y1)^2], {x1, y1, x2, y2} \[Distributed] UniformDistribution[{{0, 1}, {0, 1}, {0, 1}, {0, 1}}]]] // ComplexExpand 

2/15 + Sqrt[2]/15 + 5/48 Log[99 - 70 Sqrt[2]] - 11/192 Log[-1 + Sqrt[2]] + 173/192 Log[1 + Sqrt[2]]

confirm it. I am able to correctly calculate the integral under consideration by the change of variables which simplifies the integrand in such a way.

IntegrateChangeVariables[Inactive[Integrate][Sqrt[(x2-x1)^2+(y2-y1)^2],{x1,0,1}, {x2,0,1},{y1,0,1},{y2,0,1}],{p,q,r,s},{p==x1,q==y1,r==x2-x1,s==y2-y1}] 

Inactive[Integrate][Sqrt[r^2 + s^2], {p, 0, 1}, {q, 0, 1}, {r, -p, 1 - p}, {s, -q, 1 - q}]

Activate[ Inactive[Integrate][Sqrt[r^2 + s^2], {p, 0, 1}, {q, 0, 1}, {r, -p, 1 - p}, {s, -q, 1 - q}], Unevaluated[Integrate]] // ComplexExpand 

2/15+Sqrt[2]/15-1/64 Log[1-1/Sqrt[2]]+1/64 Log[1+1/Sqrt[2]]+13/96 Log[99-70 Sqrt[2]]-1/16 Log[-1+Sqrt[2]]+101/96 Log[1+Sqrt[2]]

Are there other ways to correctly find the integral under consideration as an iterated integral? I'd like to notice that the way to find the average length of a random interval in the unit square

in Trott, M. "The Mathematica Guidebooks Additional Material: Average Distance Distribution." http://www.mathematicaguidebooks.org/additions.shtml#S_1_14 does not consist in the calculation of an iterated integral for the integral under consideration.