## 3 Materials and Methods
### Numerical Model Implementation  
The numerical model was implemented starting from the regional-scale model developed within the MODEL-MI project, presented in the study by [PERSON] et al. (2025) [PERSON] et al. (2025). The model presented in this work was developed using the Telescopic Mesh Refinement (TMR) technique, which allows the extraction of a new model from a subregion of a larger-scale numerical model. The TMR is characterized by a rectangular domain with an extent of 4 km\({}^{2}\) (2.3 km \(\times\) 1.7 km, Figure 4), a sufficiently wide area that allows for the representation  
Figure 3: A view of the Armi pumping wells station, the MIA building’s BHE, and the hydraulic heads of the shallow aquifer from public regional data (2014). of the thermal perturbation with negligible boundary effects. Although the original model was developed in MODFLOW-2005 [57], the TMR was developed using MODFLOW-USG due to its more robust formulation, which allows the use of the CLN package and the unstructured grids (Control Volume Finite Difference formulation, CVFD). Furthermore, MODFLOW-USG allows coupled groundwater flow and heat transport simulations in the same code. The Quadtree refinement (a methodology for progressively decreasing the size of adjacent cells available for the CVFD formulation) was used in proximity to the GSHP system by defining 6 zones, starting from a maximum cell size of 16 m until a minimum cell size of 50 cm (Figure 5), defined in the area covered by the MIA building extent. The influence of the grid discretization was further analyzed by the authors in a previous study [49], showing a temperature difference lower than 0.25 \({}^{\circ}\)C 5 m downstream of the BHEs for the implemented grids. Given this result, a minimum cell size of 50 cm was assumed to be small enough for the scale and the aim of the study. The maximum grid refinement is not at the center of the quadtree rectangle because the thermal perturbation produced by BHE is expected to be downgradient of the GSHP system; therefore, being the groundwater flow directed from North to South, the need to create more refined cells towards South is essential. This technique resulted in the generation of specific refinement just in the area of interest, avoiding an unnecessary increase in the cells' number, which generally increases the computational burden. As for the horizontal discretization, a vertical spatial discretization was developed on the basis of the regional-scale model [40] (details are shown in Figure 2b). The highly detailed subsurface representation is an important advantage of this numerical model compared to those present in the literature for the Milan area (Table 1). Figure 4: A view of the numerical model domain. The implementation of 60 BHEs in MODFLOW-USG represents the first application on a real case of the adaptation of the Connected Linear Network (CLN) package, which was presented and validated in comparison with analytical solutions (Infinite Line Source and Moving Line Source) in these studies [49; 56]. The BHE is implemented as two vertical pipes (using Wells elements) and one horizontal pipe (Polyline element) to reproduce the real U-shape pipe. For the sake of simplification, the BHEs were represented as a single U, calculating the equivalent area from the actual area of the double-U BHE. The fluid inlet temperature to the BHE depends on the GSHP system operation. Here, it was assumed to be equal to 1 \({}^{\circ}\)C during the winter season and to 28 \({}^{\circ}\)C during the summer season, as reported in Table 2. \begin{table}
\begin{tabular}{c c c} \hline
**Key** & **Thickness [m]** & **Hydrogeological Unit** \\ \hline
1 & 30 & Aquifer A \\
2 & 4 & Aquifard \\
3 & 4 & Deep aquifer A \\
4 & 5 & Aquifard \\
5 & 18 & Aquifers B1 and B2 \\
6 & 5 & Aquifard \\
7 & 32 & Aquifers B3 and B4 \\
8 & 12 & Aquifard \\
9 & 40 & Aquifer C \\ \hline \end{tabular}
\end{table} 
Table 1: Vertical discretization of the numerical model. Figure 5: A view of the numerical horizontal grid discretization in the study area.