## 3 Materials and Methods
### Numerical Model Implementation
The temporal discretization of the numerical simulations was set into 4 stress periods (SPs) to reproduce the typical annual operation of a GSHP system (according to the Italian regulation [58]), for a total of three years of simulation (12 SPs). Each stress period reproduced the heating operation (during the winter season, 180 days), the cooling operation (during the summer season, 90 days), and the two intermediate seasons (spring and autumn, respectively, 50 and 40 days), as shown in Table 2. During spring and autumn, the GSHP system is not supposed to be in operation, and consequently, no flow is simulated within the BHEs in these periods. Inside each SP, the numerical solver implemented in the code is able to automatically choose the best time step size to avoid numerical dispersion into the cells where the heat transfer equation is solved. A type I boundary condition (Dirichlet condition), i.e., constant hydraulic head, was assigned along the entire perimeter of the TMR domain by imposing the hydraulic head values achieved in the last SP of the transient simulation of the calibrated regional-scale model [40]. The simulated head in the same regional model SP was used as the initial condition for the TMR. Concerning initial and boundary groundwater temperature conditions, it is generally the case within the Milan aquifers that the distribution of groundwater temperatures varies with depth. Within the Aquifer Group A (layers 1 to 4), a value of 15.2 \({}^{\circ}\)C was assigned, based on field data measured 360 m North-East of the 60 geothermal BHEs as part of a monitoring campaign performed by Politecnico di Milano. This value was imposed as a constant temperature (Dirichlet condition) at all the boundaries. As for the Aquifer Group B, a value of 14.5 \({}^{\circ}\)C was imposed, based on temperatures detected by a data logger placed inside well n. 19 of the Armi station, screened in both Aquifer Groups B and C.  
The hydrogeological properties implemented into the model were derived from the regional-scale model [40]. The hydraulic conductivity and porosity values of the aquifer implemented in the model are shown in Table 3. The specific storage value is constant in all layers and equal to 10\({}^{-4}\) m\({}^{-1}\). and the specific yield value was set equal to the porosity value. The recharge (reproducing the water infiltration on the top of the model) was assigned to the 1 st layer with values ranging between 5 \(\times\) 10\({}^{-9}\) and 4.18 \(\times\) 10\({}^{-8}\) m/s (derived from [40]) and a constant temperature value of 15.2 \({}^{\circ}\)C (assumed equal to the unperturbed temperature assumed for Aquifer A). The subsurface thermal properties were deduced by a Thermal Response Test (TRT) executed during the GSHP system design. The TRT is an in situ test performed on a single BHE aimed at evaluating the borehole thermal resistance, the ground temperature, and subsoil thermal conductivity by injecting a constant heat source into the subsurface. The thermal conductivity value estimated was equal to 2.07 W/m/K. Moreover, the maximum thermal capacity and the thermal diffusion were reported in the TRT results as 2.4 MJ/m\({}^{3}\) and 8.6 \(\times\) 10\({}^{-7}\) m\({}^{2}\)/s. The thermal dispersivity was set to zero, due to the difficulties in  
\begin{table}
\begin{tabular}{c c c c} \hline
**Layer** & **Hydraulic Conductivity [m/s]** & **Porosity [-]** & **Darcy Velocity [m/s]** \\ \hline
1 & \(2.9\times 10^{-3}\)-\(3.3\times 10^{-3}\) & 0.2 & \(7.5\times 10^{-6}\) \\
2 & \(1.8\times 10^{-6}\)-\(1.9\times 10^{-6}\) & 0.05 & \(4.5\times 10^{-9}\) \\
3 & \(3.4\times 10^{-3}\)-\(3.8\times 10^{-3}\) & 0.2 & \(8.8\times 10^{-6}\estimating this parameter for real cases. However, a value of zero allows for performing safety-oriented simulations, as the simulated groundwater temperatures may be slightly higher than the actual temperatures.