## 2 Materials and Methods
### Computational Procedure  
The governing equation, in conjunction with the boundary conditions stated earlier, is solved by means of a control-volume formulation of the finite-difference method using the open-source framework OpenFOAM [19]. A new solver, derived from the combination of the basic solvers named scalarTransportFoam and chtMultiRegionFoam, has been developed with the purpose to compute the heat transfer across adjacent solid regions with different thermophysical properties. According to the geometry of the system, fine cylindrical polar grids are used in close proximity of both the outer surface of each tube and the inner and outer surface of the insulation, whereas the remainder of the interior of the insulation and the surrounding ground at a short distance are filled with linear triangular elements. Conversely, a Cartesian grid is used to discretize the integration domain at a larger distance from the pipeline, as depicted in Figure 2.  
Starting from a uniform temperature distribution across the integration domain, the discretized governing equation is solved iteratively by way of the conjugate gradient algorithm using a second-order central difference scheme. A standard under-relaxation technique is enforced in all of the steps of the computational procedure to ensure an adequate convergence. The solution is considered to be converged when the temperature change at any grid node between two consecutive iterations is smaller than \(10^{-6}\), and the difference between the incoming and outgoing heat transfer rates, at the pipes and the ground surface, is smaller than \(10^{-2}\). Once a converged solution is achieved, the dimensionless heat transfer rates \(q_{\rm s}\) and \(q_{\rm r}\) at the outer surfaces of the supply and return pipes, normalized by the product of the temperature difference \((t_{\rm s}-t_{\rm g})\) and the insulation thermal conductivity \(k_{\rm i}\), are calculated as:  
\[q_{\rm s} =-\frac{1}{2}\int_{0}^{2\pi}\frac{\partial T}{\partial n}\Bigg{|}_ {s}d\theta \tag{7}\] \[q_{\rm r} =-\frac{1}{2}\int_{0}^{2\pi}\frac{\partial T}{\partial n}\Bigg{|}_ {r}d\theta \tag{8}\]  
On the other hand, the dimensionless heat transfer rates \(q_{\rm s}\) and \(q_{\rm r}\) can also be expressed as:  
\[q_{\rm s} =\frac{T_{\rm s}-T_{\rm g}}{\rho_{\rm pg}}+\frac{T_{\rm s}-T_{\rm r }}{\rho_{\rm pp}} \tag{9}\] \[q_{\rm r} =\frac{T_{\rm r}-T_{\rm g}}{\rho_{\rm pg}}-\frac{T_{\rm s}-T_{\rm r }}{\rho_{\rm pp}} \tag{10}\]  
where \(\rho_{\rm pg}\) and \(\rho_{\rm pp}\) are the dimensionless thermal resistances existing between each pipe and the ground surface and between the two pipes, respectively, normalized by the reciprocal of the insulation thermal conductivity \(k_{\rm i}\). Therefore, once the temperature of the return pipe is assumed to be the same as that of the supply pipe, i.e., \(T_{\rm r}=T_{\rm s}=1\), the combination of Equations (7) and (9) allows for the calculation of the value of \(\rho_{\rm pg}\). Clearly, the same result can be achieved by using Equations (8) and (10). Subsequently, once the value of \(\rho_{\rm pg}\) is replaced into Equations (9) and (10) and the temperature of the return pipe is assumed to be lower that of the supply pipe, e.g., \(T_{\rm r}=0.5\), the combination of Equations (7) and (9) allows for the calculation of the value of \(\rho_{\rm pp}\). Of course, the same result can be reached whatever value lower than unity is assumed for \(T_{\rm r}\), as well as using Equations (8) and (10) instead of Equations (7) and (9). Numerical tests on the dependence of the results on the mesh spacing, and on the extent of the whole computational domain, have been methodically performed for several combinations of the main controlling parameters. In particular, the optimal number of mesh  
Figure 2: Sketch of the discretization grid: a mesh sample of the whole integration domain (**left**) and a close-up around the interface between the insulation and the surrounding ground (**right**).