Owing to the prevailing escalation in energy costs, escalating energy demands, depletion of conventional resources, and escalating greenhouse gas emissions, there has been an expeditious upsurge in the worldwide proliferation of alternative energy sources. Geothermal systems, in particular, are highly coveted as environmentally friendly and sustainable renewable energy reservoirs, offering efficacious cooling and heating capabilities for various structures (Xu et al., 2020; Ghasemi-Fare & Basu, 2013). The energy exchange efficiency of geothermal systems is significantly contingent upon the ground temperature, enabling the discharge of energy into the ground for cooling purposes in summer and its retrieval for heating requirements during winter. However, long-term operation of geothermal energy systems can engender an imbalance in thermal recovery through the heating and cooling processes, consequently influencing the thermomechanical behavior of the subsurface over time (Baek et al., 2017). Geothermal systems, also known as ground source heat pumps (GSHPs), are designed to operate for around 20 years, while the heat source component can be utilized for an extended duration ranging from 30 to 50 years. Understanding heat transfer characteristics and long-term performance of geothermal energy systems is crucial for optimal utilization and efficiency. However, research often lacks field investigation data, specifically for long-term monitoring and validation of fluid temperatures and subsurface temperature profiles. This research conducted comprehensive long-term measurements to assess fluid and subsurface temperature behavior in a borehole heat exchanger (BHE) system. The study examined the impact of sustained BHE performance on surface and subsurface temperature dynamics. A sophisticated three-dimensional numerical model was developed based on empirical insights from the field study and validated using monitored data. The analysis included inlet/outlet fluid temperature, spatial heterogeneity of ground temperature profiles, thermal recovery processes, and alterations caused by residual thermal energy.
 Furthermore, additional temperature sensors were deployed at various depths to monitor temperature changes at distinct levels within the subground.
 In this study, an advanced three-dimensional (3D) numerical model was developed utilizing the Finite Element Method (FEM) to effectively capture the intricate heat transfer processes occurring around BHEs. The geometry mesh was generated using a specialized triangular mesh generator, enabling spatial variation of element sizes to accurately represent the geometry and temperature distribution patterns. To address the specific requirements of U-tube heat exchanger loops, where the steepest temperature gradient is anticipated, a meticulously chosen element size of 0.7 mm, was employed. As the model extends towards the lateral boundaries, particularly in the radial direction, the element size gradually increases to ensure faithful representation of the system's characteristics. To properly account for diverse vertical temperature gradients, the number of elements and the vertical distance between slices were systematically adjusted within a range of 0.01 to 5 m. This meticulous variation facilitated accurate consideration of the distinct vertical temperature profiles. Notably, the finest discretization was applied near the bottom of the BHE (approximately z = -155 m) and at the top surface, characterized by a significant temperature gradient.
 The depicted Figure 2 showcases the model's performance in analyzing the fluid profile in comparison to the monitored data over a duration of 600 hours. This comparison reveals that the maximum relative difference between the monitored and simulated values for the inlet and outlet fluid temperature does not exceed 2.649 and 2.296 degrees Celsius, respectively. These findings offer compelling evidence of the exceptional accuracy of the proposed numerical modeling approach, thus affirming its suitability for examining the long-term sustainability of BHE systems.