## 5 Conclusions  
The study demonstrates for the first time in the literature the application of the MODFLOW-USG CLN package, adapted to allow the representation in numerical models of multiple BHEs behavior in aquifers. Through the development of a 3D transient state numerical model, the subsurface thermal interference induced by a real GSHP system consisting of 60 BHEs (MIA building) was simulated and evaluated. The purpose of the study was to evaluate how the GSHP system thermal perturbation could affect the water extracted by the wells of the Armi pumping station, property of MM S.p.A., located in the city of Milan. The results of the numerical simulations highlighted that the presence of 60 BHEs belonging to the GSHP system causes moderate thermal interference in one observation point placed 80 m downgradient of the borefield (maximum 4 \({}^{\circ}\)C). Therefore, possible interferences must be considered if new GSHP systems were designed in this area, accounting for a potential GSHP system energy performance reduction. Groundwater temperature variations were more pronounced in terms of cooling the aquifer than in terms of heating, since the systems usually operate for a longer period and with greater temperature changes during the winter season. The cold thermal plume extended 2 km downstream, but with induced temperature variations lower than 1 \({}^{\circ}\)C compared to unperturbed groundwater temperature. Concerning the Armi pumping station, only limited thermal interferences (maximum of about 1.2 \({}^{\circ}\)C) were observed at the pumping well locations. Currently, such variations are expected to affect only 2 (wells n. 5 and n. 8) of the 18 wells present in the area. Therefore, no environmental or quality issues related to the influence of the GSHP system were detected for the Armi pumping station. In general, concerning the potential of shallow geothermal resources in the city of Milan, the numerical simulations highlighted some important issues related to the management of groundwater resources. Large GSHP systems such as that of MIA could determine moderate temperature variations in Aquifer Groups A, B, and C (greater than 1 \({}^{\circ}\)C), no more than 200 m along the preferential direction of groundwater flow. This distance could be taken as the first reference to protect groundwater and shallow geothermal resources. Furthermore, the presence of existing water supply infrastructures should be carefully addressed, even considering superimposition effects with existing geo-exchange systems. Therefore, in urban areas such as the city of Milan, it is becoming essential to evaluate aquifer thermal perturbations generated by geo-exchange systems using appropriate local-scale numerical models, aimed at optimizing the energy performance of GWHP or GSHP systems, avoiding thermal short-circuit effects and interactions with existing infrastructures. The following supporting information can be downloaded at [[https://www.mdpi.com/article/10.3390/app15105522/s1](https://www.mdpi.com/article/10.3390/app15105522/s1)]([https://www.mdpi.com/article/10.3390/app15105522/s1](https://www.mdpi.com/article/10.3390/app15105522/s1)), Figure S1: Thermal perturbation (plumes) induced by GSHP system in Aquifer B at the end of heating period of 1 st year simulation; Figure S2: Thermal perturbation (plumes) induced by GSHP system in Aquifer B at the end of cooling period of 1 st year simulation; Figure S3: Thermal perturbation (plumes) induced by GSHP system in Aquifer B at the end of heating period of 2 nd year simulation; Figure S4: Thermal perturbation (plumes) induced by GSHP system in Aquifer B at the end of cooling period of 2 nd year simulation; Figure S5: Thermalperturbation (plumes) induced by GSHP system in Aquifer B at the end of heating period of 3 rd year simulation; Figure S6: Thermal perturbation (plumes) induced by GSHP system in Aquifer B at the end of cooling period of 3 rd year simulation. Conceptualization, M.A., L.A., and S.B. ; methodology, M.A. and S.B. ; software, M.A. ; P.M., and S.B.