Prediction of Dynamic Temperature and Thermal Front in a Multi-Aquifer Thermal Energy Storage System with Reinjection

: The accurate temperature and thermal front prediction in aquifer thermal energy storage systems during reinjection are crucial for optimal management and sustainable utilization. In this paper, a novel two-way fully coupled thermo–hydro model was developed to investigate the dynamic thermal performance and fronts for multiple aquifer thermal energy storage systems. The model was validated using a typical model, and the evolution characteristics of wellbore temperature before and after the breakthrough of the hydraulic front and thermal front were deeply studied. Sensitivity analysis was conducted to delineate the inﬂuence of various reservoir and reinjection factors on the thermal extraction temperature (TET). The results revealed that thermal conductivity signiﬁcantly impacts the thermal extraction rate among the various reservoir factors. In contrast, volumetric heat capacity has the weakest inﬂuence and negatively correlates with the TET. Concerning the reinjection factors, the effect of the reinjection volume rate on the TET was signiﬁcantly more signiﬁcant than the reinjection temperature. Furthermore, the correlation between the TET and different properties was observed to be seriously affected by the exploitation period. The coupled model presented in this study offers insight into designing the exploitation scheme in deep reservoirs and geothermal resources.


Introduction
Geothermal energy is considered as one of the most promising sustainable, renewable, and clean energy resources, widely used for electricity generation, space heating, and greenhouse applications [1,2].In recent years, aquifer thermal energy storage (ATES) has emerged as a promising alternative for preserving excess thermal energy for future demand, owing to its advantages of enormous storage capacity, easy implementation, environmental friendliness, and good economic benefits [3][4][5].Reinjection has become a standard ATES system management strategy, which helps to enhance energy extraction efficiency and maintain reservoir pressure [6].Accurately predicting temperature behavior and front propagation in ATES reinjection systems is vital for various applications, making it crucial for decision makers and geothermal engineers in terms of optimal design and sustainable utilization.
Field experiments are costly to perform [7].However, it is crucial to determine the long-term performance of ATES systems before reinjection.Simulation technology of fluid flow and heat transfer processes has been integral to scientific research and applications over the past decades [8].Despite the many efforts made to develop accurate numerical Energies 2023, 16, 7358 2 of 24 models, even with many simulations of geothermal systems being widely applied in over 100 fields worldwide [6], the complex nature of this issue still motivates several studies [9].For example, how to accurately obtain the transient changes in various physical properties and effectively couple the seepage and temperature fields.
Significant advancements have been made in temperature modeling for fluid flow and heat transfer since the pioneering work of Lauwerier et al., who focused on simulating temperature distribution during the long-term injection of hot fluid to improve recovery percentages [10].Gringarten and Sauty et al. then presented a non-steady state mathematical model to describe transient temperature behavior during the reinjection of heat-depleted water into aquifers, building upon Lauwerier, Carslaw, and Jaeger's one-dimensional mass flow conceptual study [10][11][12][13].Chen and Reddell developed an analytical solution to describe temperature behavior during thermal injection into a confined aquifer, whether steady or unsteady [14].Their study also discussed temperature distribution in aquifer thermal energy storage (ATES) and petroleum engineering.In the 1990s, Molson et al. proposed a 3D finite element fully coupled model for simulating low-temperature groundwater flow and thermal energy transport [15].Mongelli and Pagliarulo provided a simple model of groundwater flowing parallel to the terrestrial heat flow, allowing for calculating temperature distribution and determining the extent and magnitude of the zone of influence.[16].Cheng et al. derived a multi-dimensional heat flow model for heat extraction by circulating water in a geothermal reservoir [17].Boyadjiev et al. investigated temperature profiles in a porous medium and extended the Lauwerier model [10] to include fractional extensions [18].Yang and colleagues proposed steady-state/transient temperature models that use Laplace transforms to describe ATES temperature behavior [19][20][21][22].More recently, Lee derived a coupled hydrogeological numerical model for groundwater flow and heat transport, enabling the analysis of ATES thermal behavior and recovery temperatures [1].Ganguly created a precise integral solution for a quasi-2D heat flow model in a complex geothermal reservoir.It considers advection, longitudinal conduction, and heat transfer to surrounding rock layers driven by vertical temperature differences [23].
In addition, numerical simulation software is widely used in this field.However, controlling the parameters and computational procedures can be complex, similar to a black box operation [24][25][26].Computational schemes that are faster and more accurate have been developed in reservoir simulation, while other numerical simulation codes are still under development [7].The term "black box" refers to situations in numerical simulation where a model's internal mechanisms are unknown to the user.Users can only input data and observe output results without understanding the specific operational processes within the model.In such cases, it becomes challenging for users to assess the accuracy and reliability of the model.Consequently, issues are related to low interpretability, validation difficulties, and debugging challenges.There are various types of software available for numerical simulation, including commercial software (such as ANSYS and COMSOL), open-source software (like OpenFOAM), and general-purpose software (such as MATLAB and Python), among others.Commercial software can provide various functionalities and professional technical support, making it suitable for addressing various complex problems.However, its high cost may render it less practical for small teams and researchers.On the other hand, open-source software is advantageous due to its cost-free nature, open-source code, and strong community support, offering a high degree of flexibility.Nevertheless, it often comes with a steep learning curve, necessitating users to possess a certain level of programming and numerical computing expertise.General-purpose software like MATLAB and Python offers flexibility and can be employed for diverse mathematical modeling and simulation tasks, benefiting from extensive scientific computing libraries.However, they may not perform as well as specialized finite element software in handling complex problems in specific domains.Despite these advances, few studies have been conducted on simultaneously estimating temperature profiles and predicting hydrodynamic and thermal fronts.The hydrodynamic front is where the injected fluid moves at a particular moment in ATES; the thermal Energies 2023, 16, 7358 3 of 24 front represents the position where the fluid interface temperature equals the injected fluid temperature at a specific moment.For example, in the geothermal doublet development process, one well is used for injection, while another is used for production.The reinjection of fluids tends to be much colder than the ATES system, which can cool the fluid produced from nearby wells.The premature breakthrough of hydrodynamic and thermal fronts can also cause operational problems, such as plant output running below design, added makeup wells costs, and field operations modifications [27].Therefore, it is essential to determine the propagation of hydrodynamic and thermal fronts in an ATES system.This study will establish a fully coupled model to qualitatively evaluate the relationship between temperature and hydrodynamic/thermal front.Dynamic sensitivity analysis of production temperature was carried out from the space-time perspective.
Several researchers have studied the advancement of hydrodynamic and thermal fronts in ATES through field experiments or theoretical approaches.Bodvarsson's early work in porous media highlighted two key findings: (a) the thermal front lags behind the hydrodynamic front by a constant related to volumetric heat capacities, and (b) a sudden temperature change occurs at the thermal front [28].Later, Molz et al. performed two experiments to study heat storage in aquifers and utilized test data to calibrate mathematical models that account for the interconnected processes of fluid flow and heat transfer in ATES systems [29,30].
Lippmann, Satman, Gringarten, Tsang, and others delved into the advancement of hydrodynamic and thermal fronts in various contexts, including fractured porous media.They recognized that numerical codes are essential for modeling complex, heterogeneous systems [12,[31][32][33].Researchers also examined multiphase geothermal systems with varying formation properties, developing numerous numerical codes [34][35][36].Stopa and Wojnarowski introduced a solution and equation that incorporate fluid thermal properties changing with temperature to estimate the velocity of the thermal front in geothermal reservoirs [37].Ascencio developed a theoretical temperature model for cold water reinjection in naturally fractured systems, revealing several parameters' influence on the thermal front movement rate [38].In more recent work, Abbasi and their research team proposed an exact analytical solution to describe transient temperature profiles and the advancement of the thermal front [6,39].Additionally, Shi developed a coupled thermal-hydraulic-mechanical model to investigate heat extraction performance and thermal breakthrough [2].
Although great successes have been achieved for decades, various simulation technology has been used to model ATES's transient temperature distribution and fronts analysis.Most of the mentioned studies ignore the impacts of the complex wellbore configuration and multiple adjacent formations, leading to a lack of consideration for the thermal interaction involving the wellbore, adjacent formations, and the Multiple Aquifers Thermal Energy Storage (MATES) system.The depth-dependent downward velocity of tubing fluid and the temperature-dependent thermal properties of downhole fluid were also not highly valued.In addition, few scholars have studied the relationship between temperature and fronts qualitatively or quantitatively.Thus, there is little research relating temperature inside the tubing with the MATES system's hydrodynamic front or thermal front.This study will provide insight into designing the exploitation scheme in deep reservoirs and geothermal resources.

Heat Transfer Model
Figure 1 illustrates the 2D temperature model in cylindrical coordinates for the MATES system during cold water reinjection (e.g., injecting fluids into sandstone formations with large porosity and permeability), and the heat transfer process is depicted in the schematic diagram.To facilitate the research, the complex configuration of the MATES system is divided into three layers in the horizontal direction: the tubing with fluid flow, the tubingcasing annulus with static fluid, and the tubing wall, casing, cement ring, and formations that are divided into multiple layers based on different casing programs.Moreover, in the vertical direction, the MATES system can be divided into two parts: impermeable for-mations, including the caprock, the underlying rock, and the adjacent rocks, and aquifers.The heat transfer characteristics can be summarized as: (a) The ATES region is primarily characterized by horizontal thermal convection; (b) The region inside the tubing primarily involves vertical thermal convection; (c) Other areas and boundary locations are characterized by heat conduction.
and boundary locations are characterized by heat conduction.
In addition, in this fully coupled thermal-hydraulic model, the major assumptions and approximations are as follows [8,40,41]:

Two-Way Fully Coupled Mathematical Model
Based on the heat flow process of the model, the governing equation for the temperature field was derived.This research used the implicit finite difference method to solve the transient temperature distribution, and the different discrete schemes of the model are presented in Appendix A. Furthermore, the downhole mesh generation method is also illustrated in Appendix A.

Two-Way Fully Coupled Mathematical Model
Based on the heat flow process of the model, the governing equation for the temperature field was derived.This research used the implicit finite difference method to solve the transient temperature distribution, and the different discrete schemes of the model are presented in Appendix A. Furthermore, the downhole mesh generation method is also illustrated in Appendix A.

Heat Transfer Model Inside the Tubing
During the reinjection process, the dominant heat transfers of the tubing fluid are considered to be continuous medium conduction and convection.Moreover, while the vertical component of the tubing fluid's downward seepage velocity is taken into consideration, the radial part is neglected.Therefore, the heat transfer equation inside the tubing can be expressed as [8]: where T is the temperature at a specific reinjection time, λ w is the thermal conductivity of the water, ρ w is the density of water, c w is the specific heat capacity of water, and v z denotes the downward seepage velocity of the tubing fluid.
The downward seepage velocity of the tubing fluid is proportional to volume flow; thus, the tubing fluid's downward seepage velocity (v z ) can be defined as: where q in is the total volumetric reinjection rate, r ti is the inside radius of the tubing, P A is the ratio of the volumetric reinjection rate of ATES A to the total volumetric reinjection rate, H 1 is the depth of the bottom of the ATES A, and H 2 is the depth of the bottom of the ATES B.

Heat Transfer Model from the Impermeable Formations to the Tubing Wall
For impermeable formations, heat conduction dominates heat transfer rather than convection.In addition, the tubing-casing annulus is assumed to fill with stagnant water, and thermal radiation effects are ignored.Thus, the transient temperature field governing equation for tubing, annulus, casing, cement, and impermeable formations can be given by [42]: 1 r where where λ * s is the equivalent thermal conductivity of impermeable media, (ρc) * s is the equivalent volumetric heat capacity of impermeable media, λ s is the thermal conductivity of impermeable media, λ w is the thermal conductivity of the water, ∅ is the porosity, ρ w is the density of water, c w is the specific heat capacity of water, ρ s is the density of impermeable media, and c s is the specific heat capacity of the impermeable media.

Heat Transfer Model of the MATES System
Heat conduction and convection dominated the heat transfer of the MATES system.Thus, the temperature field equation of the MATES system can be expressed as [8]: Energies 2023, 16, 7358 (ρcv) r = ρ w c w v wr (r) where λ * m is the equivalent thermal conductivity of the aquifers, (ρc) * m is the equivalent volumetric heat capacity of aquifers, λ m is the thermal conductivity of the aquifers, ∅ is the porosity, ρ w is the density of water, c w is the specific heat capacity of water, ρ m is the density of matrix in ATES, c m is the specific heat capacity of the matrix in ATES, and v wr (r) is the radial seepage velocity of the fluid in ATES.
The fluid flows of the MATES system are considered planar radial flow.Therefore, the fluid seepage velocity (v wr ) in Equation ( 9) is the function of the radial position: where v wr (r) is the radial seepage velocity of the fluid in ATES, q in is the total volumetric reinjection rate, r is the radial distance from the borehole axis, h m is the aquifer thickness, and ∅ is the porosity.

Initial and Boundary Conditions of the Model
Downhole initial (undisturbed) geothermal temperature is assumed to be a known function of depth.Thus, the initial condition can be written as: where T sur is the initial surface temperature, g T is the undisturbed geothermal gradient, and z is the vertical depth from surface.
In the radial direction, the borehole axis features an adiabatic boundary condition described as Equation ( 12); the coupled boundary conditions at interfaces between the tubing fluid and the tubing wall, the tubing fluid and the MATES system can be given by Equations ( 13) and ( 14), respectively.The outer boundary of the downhole temperature is assumed to be equal to the (undisturbed) geothermal temperature.Therefore, the outer boundary condition of the model can be defined by Equation (15).
T| r=r e = T sur + g T z where h c = 0.023 where h c is the convection heat transfer coefficient, T bf is the fluid temperature at the tubing/tubing wall interface, T bw is the tubing wall temperature at the tubing/tubing wall interface, r ti is the inside radius of the tubing, T tm is the temperature of the ATES at the interface between the tubing and the ATES, T tw is the tubing fluid temperature at the interface between the tubing and the ATES, r e is the radial distance from the borehole axis to the reservoir limit, L is the characteristic length, R ef is the Reynolds number, P rf is the Prandtl number, and n = 0.4 for the fluid being heated, n = 0.3 for the fluid being cooled.
In the longitudinal direction, we assumed that no heat transfer is occurring at the model's top and bottom; thus, their boundary conditions can be described as Equations ( 17)- (20), respectively.In addition, the coupled boundary condition at interfaces between the aquifers and adjacent formations can be defined by Equation (21).
where r ti is the inside radius of the tubing, r e is the radial distance from the borehole axis to the reservoir limit, z max is the total well depth from surface, z b is the depth of boundaries between the aquifers and the adjacent formations, T bm is the temperature of aquifers at interfaces between the aquifers and the adjacent formations, and T bs is the temperature of adjacent formations at interfaces between the aquifers and the adjacent formations.

Fluid Thermo-Physical Properties with Different Temperature
Regression analysis was used to study the thermo-physical properties of the fluid at different temperatures.The density (ρ w (T)), specific heat capacity (c w (T)), and thermal conductivity (λ w (T)) empirical formulas can be defined by Equations ( 22), ( 23) and (24), respectively.Further details on the data source and regression analysis of different parameters are presented in Appendix B.

Derivation of Hydrodynamic Front and Thermal Front Migration Model
The fluid flows of MATES system are considered as planar radial flow in this model, and the fluid seepage velocity (v wr ) of MATES system can be obtained from Equation (10).Thus, the hydrodynamic breakthrough time (t HB ) can be derived by: Based on the conservation of mass and energy theory for single-phase fluid in porous media, the thermal breakthrough time is delayed by a constant value relative to the hydrodynamic breakthrough time.Thus, it can be expressed as follows [37]: Energies 2023, 16, 7358 8 of 24 where where ∅ is the porosity, ρ m is the density of matrix in ATES, c m is the specific heat capacity of the matrix in ATES, ρ w is the density of water, and c w is the specific heat capacity of water.
Incorporating Equation ( 25), the radial position (r HF ) of the hydrodynamic front with different reinjection times can be given by: where q in is the total volumetric reinjection rate, t is the reinjection time, ∅ is the porosity, h m is the aquifer thickness, and r ti is the inside radius of the tubing.
Although any contribution from adjacent strata is neglected, the radial position (r TF ) of the thermal front at different reinjection times can be computed using the equation proposed by Ascencio et al. [38]: where where ∅ is the porosity, ρ w is the density of water, c w is the specific heat capacity of water, and (ρc) * m is the equivalent volumetric heat capacity of aquifers.

Model Validation
The proposed coupled model was programmed and compared with Li's classical analytical model [20] to ensure the mathematical equations were accurately modeled without errors.The analytical model of Li and the proposed model were based on similar concepts and reinjection forms, making the simplified model (Figure 2a) suitable for verifying the accuracy of the mathematical derivation of the proposed model.The simplified model consisted of the wellbore, an aquifer with overlying and underlying rock, and its geometry parameters, ambient temperature, and reinjection conditions are given in Table A1 (Appendix C).The thermo-physical properties of both models were derived from the same data set, provided in Table A2 (Appendix C).
Transient ATES temperature profiles of the proposed model and the analytical solution are presented in Figure 2b, with the temperature distribution trend matching well; this indicates the feasibility and accuracy of the mathematical derivation of the proposed thermal-hydraulic coupled model.The temperature difference between Li's analytical model and the proposed model can be attributed to two main reasons: (1) the proposed model of the interconnected system considered the thermal interaction between the wellbore, impermeable formations, and ATES, such as radial heat transfer of wellbore fluids and longitudinal heat transfer of ATES; however, the analytical model emphasized mathematical derivation and solution without considering the above coupled heat transfer processes.
(2) The proposed model highly valued the temperature-dependent thermo-physical properties of the fluid, while the analytical model assumed the thermo-physical properties of the fluid were spatially and temporally invariant.

Results and Discussion
Based on the discussion above, numerical solutions were utilized to simulate the temperature distribution and propagation of the hydrodynamic/thermal front.By combining the transient temperature and thermal front, the characteristics of wellbore temperature variation before and after thermal breakthrough were studied and discussed in detail.Furthermore, sensitivities were analyzed in two main categories: (a) reinjection parameters and (b) geological factors.The coefficients of variation of the outlet temperatures corresponding to various sensitivity parameters were also calculated.Appendix C includes the basic parameters used for calculation (see Tables A3-A5), and the fluid's transient thermo-physical parameters (ρw, cw, and λw) were obtained by dynamically invoking Equations ( 22)-( 24) described in Appendix B.

The Advancement Characteristic of the Hydrodynamic Front and Thermal Front
Using the model parameters of this study (Appendix C) and Equations ( 25)- (30), the radial position of the hydrodynamic and thermal fronts at different reinjection times is

Results and Discussion
Based on the discussion above, numerical solutions were utilized to simulate the temperature distribution and propagation of the hydrodynamic/thermal front.By combining the transient temperature and thermal front, the characteristics of wellbore temperature variation before and after thermal breakthrough were studied and discussed in detail.Furthermore, sensitivities were analyzed in two main categories: (a) reinjection parameters and (b) geological factors.The coefficients of variation of the outlet temperatures corresponding to various sensitivity parameters were also calculated.Appendix C includes the basic parameters used for calculation (see Tables A3-A5), and the fluid's transient thermo-physical parameters (ρ w , c w , and λ w ) were obtained by dynamically invoking Equations ( 22)-( 24) described in Appendix B.

The Advancement Characteristic of the Hydrodynamic Front and Thermal Front
Using the model parameters of this study (Appendix C) and Equations ( 25)-( 30), the radial position of the hydrodynamic and thermal fronts at different reinjection times is illustrated in Figure 3.The hydrodynamic and thermal fronts continue to advance with reinjection time.Still, the radial position of the fronts has a nonlinear relationship with reinjection time (Figure 3a).In the early injection stage, the fronts advance faster.In contrast, in the late stage of reinjection, the front migration is relatively slow.It is observed that the hydrodynamic breakthrough time and thermal breakthrough time are t HB = 406.53days and t TB = 1592.76days, respectively.Additionally, under the same conditions, we calculated the hydrodynamic front's and thermal front's breakthrough times at different reinjection volume rates (Figure 3b).The breakthrough times of the hydrodynamic and thermal front also showed a nonlinear relationship with the reinjection volume rate.Under the condition that q in < 1000 m 3 /day, the hydrodynamic and thermal front's breakthrough times vary significantly with the reinjection volume rate.When q in > 1000 m 3 /day, the change in the reinjection volume rate has little influence on thermal breakthrough.Furthermore, it is widely accepted that the location of the thermal front corresponds to the maximum radial temperature gradient.To investigate the suitability of the thermal front model derived in this study, the radial positions of fronts at different times were calculated using the above method.Figure 4 illustrates the ATES radial temperature gradient with varying times of reinjection before the hydrodynamic breakthrough.Notably, the radial position (Rmax) of the maximum value of the radial temperature gradient correlates well with the calculated thermal front position (rTF).The relative error (Derror) between rTF and Rmax can be expressed as Derror = |(Rmax − rHF)|/rHF.The relative errors for reinjection times of 5, 20, 45, and 100 days are 0.70%, 1.88%, 4.23%, and 5.94%, respectively.The detailed relative errors are presented in Figure 4b, where an increase in reinjection time leads to an increase in error.The most significant source of error may be the grid size (larger mesh sizes away from the wellbore) since the radial temperature gradient of thermal storage is based on the temperature values of grid nodes.Therefore, the front migration locations calculated using Equations ( 25)-( 30) derived in this paper are relatively accurate, enabling further research on the relationship between Furthermore, it is widely accepted that the location of the thermal front corresponds to the maximum radial temperature gradient.To investigate the suitability of the thermal front model derived in this study, the radial positions of fronts at different times were calculated using the above method.Figure 4 illustrates the ATES radial temperature gradient with varying times of reinjection before the hydrodynamic breakthrough.Notably, the radial position (R max ) of the maximum value of the radial temperature gradient correlates well with the calculated thermal front position (r TF ).The relative error (D error ) between r TF and R max can be expressed as Derror = |(Rmax − rHF)|/rHF.The relative errors for reinjection times of 5, 20, 45, and 100 days are 0.70%, 1.88%, 4.23%, and 5.94%, respectively.The detailed relative errors are presented in Figure 4b, where an increase in reinjection time leads to an increase in error.The most significant source of error may be the grid size (larger mesh sizes away from the wellbore) since the radial temperature gradient of thermal storage is based on the temperature values of grid nodes.Therefore, the front migration locations calculated using Equations ( 25)-( 30) derived in this paper are relatively accurate, enabling further research on the relationship between transient temperature in tubing and the hydrodynamic/thermal front migration position, which can provide a basis for predicting the front using wellbore temperature.(1) The temperature profile can easily determine the location of the ATES, as the spatial downhole temperature distribution resembles the shape of a "compass"; this is because the reinjection fluid of the ATES constantly migrates from the wellbore to the thermal reservoir boundary, leading to differences in cooling time for different radial positions of the ATES, resulting in varying temperatures.The temperature of the ATES near the wellbore decreases significantly, while the temperature of the ATES further from the wellbore decreases only slightly.(2) The ATES and the impermeable formations above and below are continuously cooled by the injected fluid, causing a gradual temperature decrease.As the hydrodynamic front advances, the range of cooled impermeable formations increases gradually.At the later stage of reinjection, the temperature of the ATES (1030 m ≤ z ≤ 1035 m, 1065 m ≤ z ≤ 1070 m) approaches that of the injected fluid.Moreover, Figure 6g-i reveals that during the entire reinjection life, the temperature of the ATES decreases gradually, but not indefinitely, with a "steady state" eventually reached.(3) When the hydraulic front breaks through (Figure 6f), the temperature of the heat storage boundary (R e = 518 m) does not decrease significantly.However, the boundary temperature drops significantly when the thermal front breaks through.Additionally, the entire thermal reservoir (1030 m ≤ z ≤ 1035 m, 1065 m ≤ z ≤ 1070 m) is wholly cooled for about 3000 days (Figure 6i), about twice the thermal front breakthrough time, after which, the whole heat reservoir temperature is almost equal to that of the injected fluid.
Figure 6 illustrates the evolution of the downhole temperature distribution (1000 m ≤ z ≤ 1100 m, −518.93 m ≤ r ≤ 518.93 m) during the 3000-day reinjection period.The upper part of each subgraph (Figure 6a-i) displays the temperature distribution (T), while the lower part shows the temperature variation (∆T = T − T0).The following observations can be made: (1) The temperature profile can easily determine the location of the ATES, as the spatial downhole temperature distribution resembles the shape of a "compass"; this is because the reinjection fluid of the ATES constantly migrates from the wellbore to the thermal reservoir boundary, leading to differences in cooling time for different radial positions of the ATES, resulting in varying temperatures.The temperature of the ATES near the wellbore decreases significantly, while the temperature of the ATES further from the wellbore decreases only slightly.(2) The ATES and the impermeable formations above and below are continuously cooled by the injected fluid, causing a gradual temperature decrease.As the hydrodynamic front advances, the range of cooled impermeable formations increases gradually.At the later stage of reinjection, the temperature of the ATES (1030 m ≤ z ≤ 1035 m, 1065 m ≤ z ≤ 1070 m) approaches that of the injected fluid.Moreover, Figure 6g-i reveals that during the entire reinjection life, the temperature of the ATES decreases gradually, but not indefinitely, with a "steady state" eventually reached.(3) When the hydraulic front breaks through (Figure 6f), the temperature of the heat storage boundary (Re = 518 m) does not decrease significantly.However, the boundary temperature drops significantly when the thermal front breaks through.Additionally, the entire thermal reservoir (1030 m ≤ z ≤ 1035 m, 1065 m ≤ z ≤ 1070 m) is wholly cooled for about 3000 days (Figure 6i), about twice the thermal front breakthrough time, after which, the whole heat reservoir temperature is almost equal to that of the injected fluid.

Sensitivity Analysis of ATES Thermal Extraction
To investigate the effects of different factors on thermal extraction temperature, we conducted a sensitivity analysis to provide a spatial and temporal dynamic evaluation of various parameters, such as reinjection volume rate (qin), reinjection temperature (Tin), volumetric heat capacity ((ρc)m), thermal conductivity (λm), porosity (∅m), and aquifer thickness (hm).
To compare the influence characteristics under the same conditions, we normalized Energies 2023, 16, 7358 14 of 24

Sensitivity Analysis of ATES Thermal Extraction
To investigate the effects of different factors on thermal extraction temperature, we conducted a sensitivity analysis to provide a spatial and temporal dynamic evaluation of various parameters, such as reinjection volume rate (q in ), reinjection temperature (T in ), volumetric heat capacity ((ρc) m ), thermal conductivity (λ m ), porosity (∅ m ), and aquifer thickness (h m ).
To compare the influence characteristics under the same conditions, we normalized the sensitive parameters, including q in , T in , (ρc) m , λ m , ∅ m , h m , and thermal extraction temperature (T e ).Definitions for the various normalized parameters are presented in Equations ( 31)- (37).
T enorm = T e /T a (31) where T enorm is the normalized thermal extraction temperature, T e is the thermal extraction temperature, and T a is the initial ambient temperature.
The ATES thermal extraction temperature we focused on was at a depth of 1032.6 m.To calculate the ambient temperature (T a ), we used the formula Ta = T su r + g t Z, where T sur represents the surface temperature, g t is the geothermal gradient, and Z is the depth of interest.Assuming Tsur to be 20 • C and g t to be 0.03 m/ • C, we can compute T a as Ta = 20 Figure 7 presents the normalized extraction temperatures for different normalized values of various parameters at different reinjection times.We conducted a linear regression analysis to study the relationship between normalized thermal extraction temperature and different normalized values of multiple parameters.Additionally, we defined "k slope " as the slope of the fitting curves.Our observations indicate that all parameters, including T in , (ρc) m , λ m , ∅ m , and hm, were positively correlated with thermal extraction temperature, except for qin.As reinjection time increased, the sensitivity of (ρc) m to thermal extraction temperature decreased, unlike other parameters, including T in , q in , λ m , ∅ m , and h m .Among all the sensitivity factors studied, only the sensitivity of (ρc) m to thermal extraction temperature decreased with reinjection time.K slope fell from 0.022 to 0.0155 when reinjection time varied from 1.0 to 10.0 years, indicating that the impact of (ρc) m on thermal extraction temperature can be disregarded in the later stages of reinjection time.Additionally, the K slope of q in was significantly greater than that of T in for reinjection factors.
(ρc)m to thermal extraction temperature decreased, unlike other parameters, including Tin, qin, λm, ∅m, and hm.Among all the sensitivity factors studied, only the sensitivity of (ρc)m to thermal extraction temperature decreased with reinjection time.Kslope fell from 0.022 to 0.0155 when reinjection time varied from 1.0 to 10.0 years, indicating that the impact of (ρc)m on thermal extraction temperature can be disregarded in the later stages of reinjection time.Additionally, the Kslope of qin was significantly greater than that of Tin for reinjection factors.

Comprehensive Comparison of Sensitivity Analysis
To comprehensively compare the influence degree of various parameters on the thermal extraction temperature in this study, based on the principle of mathematical

Comprehensive Comparison of Sensitivity Analysis
To comprehensively compare the influence degree of various parameters on the thermal extraction temperature in this study, based on the principle of mathematical statistics, we adopted the coefficient of variation as the index of sensitivity evaluation [43,44]: Energies 2023, 16, 7358 where C v is the coefficient of variation, s is the standard deviation of the sample data, |y| is the average of sample data, y i represents a data sample, and N represents the number of sample data.The coefficient of variation (C v ) represents the sensitivity of thermal extraction temperature to different influencing factors.A higher C v value indicates a more significant influence of that factor on thermal extraction temperature.Figure 8 summarizes the results of the sensitivity comparison of the parameters studied.Among injection parameters, reinjection volume rate (q in ) influences the thermal extraction temperature the most.Regarding ATES attributes, thermal conductivity (λ m ) has the highest sensitivity within the studied parameter ranges, followed by porosity (∅ m ).The reservoir volume heat capacity ((ρc) m ) and reservoir thickness (h m ) have relatively little impact on thermal extraction temperature.Therefore, thermal conductivity and porosity are essential for evaluating thermal reservoirs and extraction.In addition, the reinjection volume rate is a critical evaluation parameter among injection parameters.As reinjection time increases, the influence of different parameters on heat extraction temperature is more significant.Thus, the reinjection times significantly influence the correlations between thermal extraction temperature and thermo-physical properties.Therefore, reinjection time is also a critical factor to consider in heat storage development.
statistics, we adopted the coefficient of variation as the index of sensitivity evaluation [43,44]: (39) where Cv is the coefficient of variation, s is the standard deviation of the sample data, | ̅| is the average of sample data, yi represents a data sample, and N represents the number of sample data.
The coefficient of variation (Cv) represents the sensitivity of thermal extraction temperature to different influencing factors.A higher Cv value indicates a more significant influence of that factor on thermal extraction temperature.Figure 8 summarizes the results of the sensitivity comparison of the parameters studied.Among injection parameters, reinjection volume rate (qin) influences the thermal extraction temperature the most.Regarding ATES attributes, thermal conductivity (λm) has the highest sensitivity within the studied parameter ranges, followed by porosity (∅m).The reservoir volume heat capacity ((ρc)m) and reservoir thickness (hm) have relatively little impact on thermal extraction temperature.Therefore, thermal conductivity and porosity are essential for evaluating thermal reservoirs and extraction.In addition, the reinjection volume rate is a critical evaluation parameter among injection parameters.As reinjection time increases, the influence of different parameters on heat extraction temperature is more significant.Thus, the reinjection times significantly influence the correlations between thermal extraction temperature and thermo-physical properties.Therefore, reinjection time is also a critical factor to consider in heat storage development.

Conclusions
This study established a fully coupled model to qualitatively evaluate the relationship between temperature and hydrodynamic/thermal front.Dynamic sensitivity analysis of production temperature was carried out from the space-time perspective.The following conclusions have been drawn.

Conclusions
This study established a fully coupled model to qualitatively evaluate the relationship between temperature and hydrodynamic/thermal front.Dynamic sensitivity analysis of production temperature was carried out from the space-time perspective.The following conclusions have been drawn.
(a) The accuracy of the thermal front migration calculation model for geothermal reservoirs with radial fluid seepage has been confirmed.The position of the hydraulic front can be used to accurately predict the position of the thermal front before the hydraulic breakthrough, as the radial position (R max ) of the maximum radial temperature gradient corresponds well with the calculated radial position (r HF ) of the thermal front.The relative error between R max and r HF ranges from 0.70% to 5.94%, providing a reliable basis for forecasting the thermal front position various times before hydraulic breakthrough.
The temperature increase rate in the upper parts of ATES A (1010-1020 m) and ATES B (1050-1060 m) is almost negligible.However, the temperature increase rate turns negative after the thermal breakthrough (t > t TB = 1592.76days).Furthermore, the temperature in the bottom hole (z = 1100 m), located approximately 40 m from the base of ATES B, remains close to the formation's original temperature before the thermal breakthrough.Following the thermal breakthrough, the temperature drops by approximately 13%.These trends can serve as an indicator to detect hydrodynamic or thermal breakthroughs based on temperature monitoring inside the tubing.(c) The temperature gradient in the ATES section is close to zero (grad(T) = 0).Conversely, the temperature gradient in the impermeable formations section is more significant, allowing for identifying the location of the ATES or impermeable formations through analysis of the temperature inside the tubing.(d) The sensitivity analysis revealed that the normalized thermal extraction temperature variation with various normalized parameters is nearly linear.Furthermore, the following significant findings can be noted: (i) Among the reservoir factors (such as λ m , (ρc) m , ∅ m , and h m of the ATES), λ m holds the most significant influence on the thermal extraction rate.In contrast, (ρc) m has the weakest impact.In contrast, q in had a significantly higher impact on thermal extraction temperature than T in , with only q in showing a negative correlation with the thermal extraction temperature.(ii) Thermal conductivity and porosity are vital evaluation indices for assessing thermal reservoir development.Similarly, the reinjection volume rate is critical among the various injection parameters.As the reinjection time increases, the impact of different parameters on heat extraction temperature becomes more significant.Hence, the reinjection time must be meticulously considered in evaluating heat storage or heat extraction.(e) This study is applicable not only for single ATES but also for MATES, allowing us to consider the aquifer and adjacent formation as having multiple layers with dynamic thermo-physical properties.Moreover, the proposed coupled thermal-hydraulic model, as opposed to analytical or semi-analytical models, offers the capability to predict the complete dynamic temperature distribution within the wellbore of complex configuration and reservoirs under more realistic heat exchange conditions.(f) To make the model more closely resemble the natural geothermal geological environment, there are two main tasks for the next step: first, considering random fractures' impact on temperature and the thermal front, and secondly, studying the relationship between transient temperatures at different positions within the wellbore and the fluid front/thermal front.The ratio of the volumetric reinjection rate of ATES A to the total volumetric reinjection rate (q in ), % P rf Prandtl number, dimensionless q in Total volumetric reinjection rate, m 3 /s q norm Normalized total volumetric reinjection rate, dimensionless r The radial distance from the borehole axis, m r ai Inside radius of cement A, m r ci Inside the radius of the casing, m r e The radial distance from the borehole axis to the reservoir limit, m r HF The radial position of the hydrodynamic front with different reinjection times, m r TF The The temperature of aquifers at interfaces between the aquifers and the adjacent formations, • C T bs The temperature of adjacent formations at interfaces between the aquifers and the adjacent formations,

( a )
Undisturbed geothermal temperature of the downhole is a known function of depth, and the vertical borehole radius remains unchanged with depth.(b) The reinjection well thoroughly penetrates the MATES system, and the fluid flow is steady state.(c) MATES fluid flows are considered planar radial flow.(d) Fitting functions of the thermo-physical parameters with different temperatures are suitable for the coupled thermal-hydraulic model.

Figure 1 .
Figure 1.The schematic of the thermal-hydraulic coupled model.

Figure 1 .
Figure 1.The schematic of the thermal-hydraulic coupled model.In addition, in this fully coupled thermal-hydraulic model, the major assumptions and approximations are as follows [8,40,41]: (a) Undisturbed geothermal temperature of the downhole is a known function of depth, and the vertical borehole radius remains unchanged with depth.(b) The reinjection well thoroughly penetrates the MATES system, and the fluid flow is steady state.(c) MATES fluid flows are considered planar radial flow.(d) Fitting functions of the thermo-physical parameters with different temperatures are suitable for the coupled thermal-hydraulic model.

Figure 2 .
Figure 2. (a) The schematic representation of the simplified model.(b) Comparison of the transient ATES temperature profile of the proposed model with an analytical solution of Li's classical analytical model [20].

Figure 2 .
Figure 2. (a) The schematic representation of the simplified model.(b) Comparison of the transient ATES temperature profile of the proposed model with an analytical solution of Li's classical analytical model [20].

Energies 2023 ,Figure 3 .
Figure 3. (a) The radial position of hydrodynamic and thermal fronts during the reinjection.(b) The breakthrough times of the hydrodynamic and thermal front at varying reinjection volume rates.

Figure 3 .
Figure 3. (a) The radial position of hydrodynamic and thermal fronts during the reinjection.(b) The breakthrough times of the hydrodynamic and thermal front at varying reinjection volume rates.

Energies 2023 ,
16,  x FOR PEER REVIEW 11 of 24 transient temperature in tubing and the hydrodynamic/thermal front migration position, which can provide a basis for predicting the front using wellbore temperature.

Figure 4 .
Figure 4. ATES radial temperature gradient with different reinjection times before the hydrodynamic breakthrough time.(a-d) are the ATES radial temperature gradient profiles of 5, 20, 45, and 100 days of reinjection time, respectively.(e) The relative error between the rHF and Rmax before the hydrodynamic breakthrough time.

3. 2 .Figure 5 ( 1 )Figure 4 .
Figure 5 displays the temperature distribution inside the tubing at the tubing/tubing wall interface (r = 0.04 m, 1000 m ≤ z ≤ 1100 m) with different reinjection times.The following observations can be made: (1) The tubing fluid in sections of adjacent formations A (1000 m ≤ z ≤ 1030 m) and B (1035 m ≤ z ≤ 1065 m) is continuously cooled by the injected fluid, with heat

3. 2 .Figure 5 ( 1 )
Figure 5 displays the temperature distribution inside the tubing at the tubing/tubing wall interface (r = 0.04 m, 1000 m ≤ z ≤ 1100 m) with different reinjection times.The following observations can be made: (1) The tubing fluid in sections of adjacent formations A (1000 m ≤ z ≤ 1030 m) and B (1035 m ≤ z ≤ 1065 m) is continuously cooled by the injected fluid, with heat convection and conduction of wellbore fluid dominating the heat transfer, causing a gradual temperature decrease until it approaches the injected fluid temperature.However, for wellbore fluid in the underlying rock depth range (1070 m ≤ z ≤ 1100 m), longitudinal heat conduction dominates the heat transfer, resulting in a more gradual temperature decrease.(2) In sections of ATES A (1030 m ≤ z ≤ 1035 m) and B (1065 m ≤ z ≤ 1070 m), the fluid temperature in the tubing is equal to the reinjection temperature (Tin = 20 °C) throughout the injection life.As a result, the longitudinal temperature gradient inside

Figure 5 .
Figure 5. Temperature distribution inside the tubing at the tubing-tubing wall interface (r = 0.04 m, 1000 m ≤ z ≤ 1100 m) with different injection times.

Figure 5 .
Figure 5. Temperature distribution inside the tubing at the tubing-tubing wall interface (r = 0.04 m, 1000 m ≤ z ≤ 1100 m) with different injection times.

3. 3 .
Dynamic Thermal Performance Analysis of MATES 3.3.1.Spatial Distribution of the Downhole Transient Temperature (1000 m ≤ z ≤ 1100 m, −518.93 m ≤ r ≤ 518.93 m) Figure 6 illustrates the evolution of the downhole temperature distribution (1000 m ≤ z ≤ 1100 m, −518.93 m ≤ r ≤ 518.93 m) during the 3000-day reinjection period.The upper part of each subgraph (Figure 6a-i) displays the temperature distribution (T), while the lower part shows the temperature variation (∆T = T − T 0 ).The following observations can be made:

24 Figure 6 .
Figure 6.Evolution of the spatial distribution of the downhole temperature (1000 m ≤ z ≤ 1100 m, −518.93 m ≤ r ≤ 518.93 m) during the 1500 days.(a-i) represent the temperature distribution from 5.0 days to 3000 days.

Figure 6 .
Figure 6.Evolution of the spatial distribution of the downhole temperature (1000 m ≤ z ≤ 1100 m, −518.93 m ≤ r ≤ 518.93 m) during the 1500 days.(a-i) represent the temperature distribution from 5.0 days to 3000 days.

Figure 7 .
Figure 7.The normalized extraction temperatures for different normalized values of various parameters with varying times of reinjection.

Figure 7 .
Figure 7.The normalized extraction temperatures for different normalized values of various parameters with varying times of reinjection.

Figure 8 .
Figure 8. Coefficient of variation of the outlet temperature corresponding to the various sensitivity parameters.

Figure 8 .
Figure 8. Coefficient of variation of the outlet temperature corresponding to the various sensitivity parameters.
c mSpecific heat capacity of the matrix in ATES, J/(kg•K) c s Specific heat capacity of the impermeable media, J/(kg•K) The depth of the bottom of the ATES A, m H 2The depth of the bottom of the ATES B, m nThe coefficients in Equation (16), n = 0.3 K slopeThe slope of normalized extraction temperature with different normalized values of various parameters, dimensionless P A radial position of the thermal front with different reinjection times, m

Table A2 .
Data used in model comparison.
• C

Table A4 .
Thermo-physical properties of various media.

Table A5 .
Ambient temperature and reinjection conditions.