Quantitative analysis of the influence of middle-deep structural characteristics on shallow geothermal field

The study of the shallow geothermal field can provide precise parameters for the optimization design of ground source heat pump systems, and offer scientific basis for urban energy planning and the selection of new energy stations in new areas. The middle and deep structure has a certain influence on the shallow geothermal energy, and the influence law of different structures (bedrock characteristics, fault structure) is quite different. Based on the geological structure and stratigraphic characteristics of Gaoyang geothermal field in Xiong 'an New Area, the theoretical geological model and mathematical model of different structural characteristics are constructed. The influence of bedrock burial depth and fault cutting depth on shallow geothermal field is quantitatively analyzed. The study shows that the sensitivity of the fault structure to the temperature field is about 1.9 times that of the uplift structure in the simulation of the heat flow coupling field. The contribution rate coefficient increases with the increase of temperature field depth, and the contribution rate of the fault structure to the shallow temperature field in the heat-flow coupling mode is about 13 times that of the thermal-solid coupling. The shallower the buried depth of the uplift is, the more intense its influence on the shallow temperature field will be. The influence of fault structures on the temperature field is much stronger than that of uplift structures. The seepage field of the uplift structure and the fluid in the fault structure directly affect the shallow geothermal field.

Introduction

Reducing carbon dioxide emissions has become a global concern. Shallow geothermal energy is a new renewable and environmentally friendly energy source with low cost, wide coverage, low pollution and great development potential (Gude, 2016; Urchueguia, 2016; Pang et al., 2018; Roka et al., 2023). Large-scale development and utilization of shallow geothermal energy is the most effective measure to alleviate the energy crisis and achieve national heating and cooling decarbonization (Lund and Toth, 2021; Zhang, 2022; Walch et al., 2022; Ahmed et al., 2022; Figueira et al., 2024). Due to the significant differences in land use, underground space management, climate, geography, geology and hydrological conditions in each city, the utilization of shallow geothermal energy varies from city to city (Xu et al., 2020; Wang G. L. et al., 2024). Studying the characteristics of the shallow geothermal field can guide the optimization application of ground source heat pump technology and assess the potential of regional shallow geothermal resources (Li et al., 2017; Luan et al., 2013).

The regional distribution of shallow geothermal energy is influenced by multiple factors, including the characteristics of bedrock, geothermal gradient, lithologic structure, features of fracture structure, and the distribution characteristics of groundwater (Wei et al., 2012; Wei et al., 2020; Wei and Li, 2021; Li et al., 2017; Zhu et al., 2022). The top boundary of bedrock and its geothermal field are the fundamental conditions that determine the spatial distribution of shallow geothermal energy. The buried depth of bedrock is positively correlated with the buried depth of the constant temperature layer roof. Active faults are the key factors affecting and controlling shallow geothermal energy (Wei et al., 2020). From the perspective of geothermal heat transfer, active faults change bedrock fissures, thus changing the channels of fluid and heat transfer from the deep earth to the surface. Therefore, it also affects the distribution characteristics of shallow geothermal energy. The regional geological structure and the deep crustal structure exert a controlling influence on the ground temperature, determining the overall trend of the ground temperature distribution. In addition, factors such as rock properties, groundwater activity, topography and precipitation also significantly influence the specific distribution of local ground temperatures (Lei et al., 2017; Li, 2021; Yuan et al., 2024; Li et al., 2024; Lei et al., 2019). There are few studies on the influence of middle-deep structure on shallow geothermal field (Wei et al., 2012; Wei et al., 2020; Wang et al., 2017). The influence of middle-deep geological structure on shallow geothermal field is the result of the interaction and mutual influence of multiple factors. However, the main controlling factors of the influence of different geological structures including bedrock morphology, bedrock burial depth and fault morphology, cutting depth and direction on shallow geothermal field are not clear.

Numerical simulation is a very important research method. Its core role is to reproduce and expand the physical similarity simulation experiment, and to realize the quantitative analysis of key parameters that are difficult to accurately capture in the experiment, such as stress field distribution, seepage field evolution and temperature field change (Sun et al., 2025; Wang Z. T. et al., 2025; Teng et al., 2025). Geological processes are the main controlling factors for the development of geothermal systems, and numerical simulation is a fundamental method for reproducing and evaluating the processes of geothermal systems (Franco and Vaccaro, 2014; Torresan et al., 2022; José et al., 2022; Ioan and Alexandru, 2024). Many scholars have carried out numerical simulation research on the characteristics of geothermal field through numerical simulation methods. The tectonic structures are fundamental connections between the recharge area and the exploitation field (Torresan et al., 2022; Wang L. et al., 2025). The temperature field distribution varied in different strata, which was closely related to the thermal properties of the strata and the groundwater flow (Lei et al., 2019; Yue et al., 2017; Chen et al., 2021; Wang W. L. et al., 2024; Feng et al., 2024). The temperature field and water density field are strongly controlled by the fault system, and there are high temperature and low density zones within the fault system (Zeng et al., 2020). The shallow geothermal field of the crust is mainly controlled by the undulation of the basement and the thickness of the sedimentary strata (Shi et al., 2022; Darius et al., 2011). Tang (2019) conducted a numerical simulation study and found that as the depth increases, the ground temperature is mainly affected by the terrestrial heat flow. The simulation of these geothermal resources is mainly focuses on deep geothermal resources, or on the temperature field variation characteristics during the geothermal extraction process. There are relatively few numerical simulations on the internal mechanism of heat exchange between deep and shallow geothermal energy, especially the comparative research and analysis based on the thermal-fluid-solid coupling mathematical model.

There are few studies on the influence of mid-deep structure on the temperature field characteristics of shallow geothermal energy and the heat transfer mechanism between them (Wu and Song, 2010; Wang et al., 2019). Starting from the characteristics of middle and deep structures, a three-dimensional geological model with different main control structural characteristics is constructed by using the existing geophysical data, drilling data and other field data. Fully considering the energy exchange relationships among temperature, fluid and solid, a thermal-fluid-solid coupling mathematical model is constructed, and the three-dimensional data volume of complex temperature field with different main control structures is obtained by numerical simulation. The heat transfer mechanism between the deep and shallow layers in the study area is analyzed based on the two heat exchange modes of heat conduction (thermo-solid coupling) and heat convection (heat-flow coupling), and the dominant heat exchange mode in the study area is given.

In the actual stratigraphic research, it is extremely difficult to analyze and study a specific influencing factor independently, and it is impossible to carry out comparative research. In order to further study the influence of different geological structures on shallow geothermal field, the main controlling factors of the influence of middle and deep structures on shallow geothermal field are analyzed. In this study, several ideal geological models are constructed, including reference model, bedrock structure model and fault structure model. By using the numerical simulation method, the influence of bedrock buried depth and fault structure on shallow geothermal field is analyzed from the quantitative perspective.

Model construction

Conceptual model construction

The influencing factors of shallow temperature field are complex and difficult to be analyzed separately. The theoretical model can be used to analyze the influence law of different influencing factors from a single factor. The geological models in the study are divided into three groups, and the overall size of the three types of models is 10 km × 8 km × 4 km. Referring to the stratigraphic characteristics of Xiong 'an New Area, it is divided into five layers, including Pingyuan Formation, Minghuazhen Formation, Dongying Formation, Cambrian and Wumishan Formation. The first group is the reference model. The model stratum is isotropic and does not contain any other geological structure. The reference model is mainly used for comparative reference when analyzing subsequent models. The second group is the uplift structure model, and the buried depth of the uplift structure of each model is different. The third group is the fault structure model, and the fault cutting depth of each model is different (Figure 1). The design of the uplift structure model is based on the main uplift structure characteristics of Xiong 'an New Area. Two uplift structure models are designed. The buried depth of the two uplift structures is different, and the other parameters are consistent. By changing the burial depth of the uplift structures, the general law of the influence of the burial depth of uplift structures on the shallow geothermal field was studied. The buried depth of the uplift structure model 1 is set to 2,100 m, and the buried depth of the uplift structure model 2 is set to 1,400 m. The design of the fault structure model is based on the main fault structure characteristics of Xiong 'an New Area. Three types of fault structure models are designed. The cutting depth of the three fault structures is different, and the other parameters are consistent. By changing the cutting depth of the fault structure, the general law of the influence of the fault structure characteristics on the shallow geothermal field is studied. The cutting depth of fault structure model 1 is set to 2,450 m, the cutting depth of fault structure model 2 is set to 3,400 m, and the cutting depth of fault structure model 3 is set to 3,800 m. The dip angle of the fault structure of the three models is set to 45°.The cutting depth and dip angle in the model are determined based on the actual buried depth of the fault structure in Xiong 'an New Area (Wang et al., 2020; Wang G. L. et al., 2024; Ma et al., 2020; Gao et al., 2023; Dai et al., 2023).

Figure 1

Figure 1. Geological model schematic diagram (a) reference model; (b) uplift structure model (with burial depth of 2.1 km); (c) uplift structure model (with burial depth of 1.4 km); (d) fault structure model (with cutting depth of 2.45 km); (e) fault structure model (with cutting depth of 3.4 km); (f) fault structure model (with cutting depth of 3.8 km).

All the information of the geological model needs to be obtained through field exploration, and the original data of the formation can be obtained directly through logging, geochemical and geophysical exploration. Accurate and abundant original data is the basis of geological model construction. The selection and setting of the parameters of the theoretical model, the uplift structure model and the fault structure model in this study are also based on the actual stratigraphic characteristics in the field. The geological structures and the stratigraphic lithology parameters of the two models are based on the Xiong 'an New Area.

Mathematical model

When considering transient conditions and groundwater flow, the energy conservation equation used is (COMSOL Inc., 2020):

$\underset{①}{\underbrace{{\left(\rho c_p\right)}_{eff}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}}}+\underset{②}{\underbrace{\rho c_p\mathit{u}·\nabla T}}-\underset{③}{\underbrace{\nabla ·{\lambda }_{eff}\nabla T}}=Q(1)$

If there is a fault zone, the energy conservation equation for the fault zone is:

$\underset{①}{\underbrace{d_f{\left(\rho c_p\right)}_{eff}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}}}+\underset{②}{\underbrace{d_f\rho c_p\mathit{u}·{\nabla }_{T}T}}-\underset{③}{\underbrace{{\nabla }_T·\left(d_f{\lambda }_{eff}{\nabla }_{T}T\right)}}=d_{f}Q(2)$

Whereρis density (kg/m³); c_p is specific heat capacity (J/(kg·K)); T is temperature(K); t is time (s); $\mathit{u}$ is the fluid velocity (m/s); ${\lambda }_{eff}$ is effective mean conductivity (W/(m·K)); d_f is width of the fault zone(m); Q is fluid mass (kg/(m³·s)).

Flow field control equation:

Based on the theory of porous elasticity, the mass conservation equation for the fractured rock mass in geothermal reservoir (The terms ② in Equation 2):

$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\varphi {\rho }_w\right)+\nabla ·\left({\rho }_w\mathit{u}\right)=Q(3)$

where $u$ is the fluid velocity (m/s); φ is Porosity; ${\rho }_w$ is water density (kg/m³). Assuming that the Darcy's law is satisfied, then u in Equation 3 is:

$\mathit{u}=-\frac{k}{{\mu }_w}\left(\nabla p+{\rho }_w\mathit{g}\right)(4)$

where k is permeability (m²); ${\mu }_w$ is the dynamic viscosity of water (Pa·s); p is fluid pressure (Pa); g is the acceleration of gravity (m/s²).

If there is a fault zone, the fluid flow equation within the fault zone is the Equations 2, 3 have been transformed into Equation 5:

$d_f\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\varphi }_f{\rho }_w\right)+{\nabla }_T·\left(d_f{\rho }_w\mathit{u}\right)=d_{f}Q(5)$

$\mathit{u}=-\frac{k_f}{{\mu }_w}\left({\nabla }_{T}p+{\rho }_w\mathit{g}\right)(6)$

Where ${\varnothing }_f$ is porosity of the fault zone(%); $k_f$ is the fault permeability, which satisfies the cubic law and is expressed as:

$k_f=\frac{d_f^2}{12f_f}(7)$

where $f_f$ is the fracture roughness coefficient. The Equations 5–7 is the fluid flow equation within the fault zone.

Thermal-solid coupling situation:

When considering the solid deformation of the rock layer but not the groundwater flow, terms ① and ② in Equation 1 are deleted, and the energy conservation equation is simplified to:

$-\nabla ·{\lambda }_{eff}\nabla T+K{\alpha }_{T}T\frac{\mathrm{\partial }{\epsilon }_v}{\mathrm{\partial }t}=Q(8)$

If there is a fault zone, the energy conservation equation for the fault zone is:

$-{\nabla }_T·\left(d_f{\lambda }_{eff}{\nabla }_{T}T\right)+d_{f}K{\alpha }_{T}T\frac{\mathrm{\partial }{\epsilon }_v}{\mathrm{\partial }t}=d_{f}Q(9)$

Where K is bulk Modulus (Mpa); ${\alpha }_T$ is thermal Expansion Coefficient (K⁻¹); ${\mathrm{\epsilon }}_{\mathrm{v}}$ is volumetric Strain;

Solid deformation control equation

Based on the assumption of linear elasticity and small deformation, the deformation equilibrium equation of rock matrix under pure gravity is obtained as follows:

$\nabla ·\mathit{\sigma }+{\rho }_{av}\mathit{g}=0(10)$

Based on the linear elastic constitutive, the elastic stress tensor ($\mathit{\sigma }$) in Equation 10 is:

$\mathit{\sigma }=\mathbf{D}:{\mathit{\epsilon }}_{el}(11)$

When considering the thermal expansion, the elastic strain (in Equation 11) is equal to the total strain tensor minus the thermal strain tensor ${\mathit{\epsilon }}_{th}$

${\mathit{\epsilon }}_{el}=\frac{1}{2}\left(\nabla \mathit{u}+{\nabla }^T\mathit{u}\right)=\mathit{\epsilon }-{\alpha }_T\left(T-T_{ref}\right)(12)$

Substituting Equations 8, 9, 11, 12 into Equation 10, the final transformed equation can be expressed as:

$\left(\lambda +G\right)\nabla {\epsilon }_v+G{\nabla }^2\mathit{u}-K{\alpha }_T\nabla T+{\rho }_{av}\mathit{g}=0(13)$

Where: $\mathbf{D}$ is the fourth-order elastic matrix tensor, $\lambda$ and $G$ are Lamé constants, ${\epsilon }_v$ is the volumetric strain, $\mathit{u}$ is the displacement vector, ${\rho }_{av}=\varphi {\rho }_{water}+\left(1-\varphi \right){\rho }_s$ represents the average density, $T_{ref}$ is the reference temperature (K).

Numerical model

Boundary and inital conditions

According to the meteorological and geothermal geological data for many years, the depth of the constant temperature zone is determined to be 25 m. The top interface of the model is set as a constant temperature boundary, and the temperature is set to be 10 °C. The bottom interface is also set as a the constant temperature boundary, and the temperature is set to be 120 °C (Gao et al., 2023). The boundary conditions on both sides are open boundaries. In the boundary setting of seepage field, the left side of the model is set as no flow, the right side is set as a flow boundary, the bottom interface is also set as a flow boundary, and the bottom interface and the right side flow boundary are both set to 1 × 10⁻⁸ kg/m²·s. The pressure boundary at the top of the pressure boundary model is set to 0.1 MPa.

Model parameters

The geological parameters in the model are based on the field logging, drilling and rock samples of Xiong 'an New Area, ensuring that the geological parameters in the model are close to the actual geological conditions. This enables the simulation results to better reflect the actual situation and have certain guiding significance (Ma et al., 2020). The parameters of each layer of the model are shown in Table 1.

Table 1

Table 1. Model parameter table.

Model discretization

Considering the balance between calculation accuracy and computational efficiency, the boundary and fracture surface meshes are refined, and the grid division of each model is shown in Figure 2.

Figure 2

Figure 2. Model grid division diagram (a) Reference model; (b) Uplift structure model (buried depth 2.1 km); (c) Uplift structure model (buried depth 1.4 km); (d) Fault structure model (cutting depth 2.45 km); (e) Fault structure model (cutting depth 3.4 km); (f) Fault structure model (cutting depth 3.8 km).

Simulation results

Simulation results of different structures

The simulation results of different models are obtained by simulation (Figure 3) based on Equation 13. Compared with the reference model, the temperature of the uplift model 1 increases slightly at the uplift position (Figure 3b), while the temperature increase in the uplift model 2 is more obvious (Figure 3c). In the center of the uplift, the temperature reaches 85 °C at a depth of 1,500 m; By comparing the three fault models, it can be seen that the temperature at the fault location changes significantly. The heat flow at the bottom rapidly transfers to the upper part through the fault, causing the temperature to rise. The deeper the fault cutting, the higher the temperature. At the buried depth of 2 km, the temperature of fault model 1 is 80 °C (Figure 3d), the temperature of fault model 2 is 100° (Figure 3e), and the temperature of fault model 3 increases to 110° (Figure 3f). With the increase of the buried depth of the fault, the heat increases along the fault direction.

Figure 3

Figure 3. Simulation results of heat-fluid coupling (a) Temperature field of the reference model; (b) Temperature field of the uplift structure (buried depth of 2.1 km); (c) Temperature field of the uplift structure (buried depth of 1.4 km); (d) Temperature field of the fault structure (cutting depth 2.45 km); (e) Temperature field of the fault structure (cutting depth 3.4 km); (f) Temperature field of the fault structure (cutting depth 3.8 km).

Verification analysis of simulation results reliability

Although the models studied in this project are all ideal models based on rational conditions, all the model construction parameters and geological layer parameters in the models are constructed and set according to the actual field measurement data. By comparing the temperature data of the simulation results with the actual logging temperature data, the rationality of the model and parameters, as well as their similarity to the actual situation, can be verified to a certain extent.

Using the numerical simulation results obtained by simulation, the temperature curves of different models are extracted and compared with the temperature curves of wells D17, D18, D32, D34 and D35 in the Xiong 'an New Area (Figure 4). The results show that there is a certain difference between the simulated uplift structure and the actual logging temperature curve, and the overall temperature value is lower than that of the actual logging curve. This is because in the actual on-site conditions, uplift structures are often accompanied by the fault structure. Therefore, during the heat transfer process in the middle and deep layers, the heat transfer of the uplift structures in the simulation is less than that in the actual situation where both uplift structures and fault structures are present. However, the temperature change trends between the two remain quite consistent. The simulated temperature variations of the faults show good consistency with the actual logging temperature curves in terms of temperature values and temperature gradient changes. By comparing the simulation results with the actual logging temperature curves, it can be determined that the setting of the model parameters is reasonable, and the simulation results have certain reliability.

Figure 4

Figure 4. Comparison of simulated data and field logging temperature curve. The heat transfer of the uplift structures in the simulation is less than that in the actual situation. The simulated temperature variations of the faults show good consistency with the actual logging temperature curves.

Discussion and analysis

Comparative analysis of sensitivity and contribution

When analyzing the influence of different middle and deep structures on the shallow temperature field, the contribution rate parameter (M) and the sensitivity coefficient (E) are mainly analyzed and studied. The contribution rate parameters are used to study and analyze the contribution of different structures to shallow geothermal energy, and the sensitivity coefficient is used to study and analyze the sensitivity of shallow geothermal energy to different structures.

The calculation formulas for the contribution rate parameter (M) and the sensitivity coefficient (E) are as follows:

$\mathrm{M}=\Delta \mathrm{t}/{\mathrm{T}}_0(14)$

$\mathrm{E}=\left(\Delta \mathrm{t}/{\mathrm{T}}_0\right)/\left(\Delta \mathrm{h}/{\mathrm{H}}_0\right)(15)$

Where T₀ is the temperature of the initial temperature field, Δt is the shallow temperature change after adding or changing the uplift or fault characteristics; H₀ is the parameter value (depth) of the uplift or fault in the initial state, and Δh is the change of the parameter.

The temperature curves of uplift structure and fault structure model at different buried depths are extracted (Figure 5; Table 2). The sensitivity coefficient and contribution rate parameters of the uplift structure and the fault structure are calculated by using the extracted temperature curve data (Table 3).

Figure 5

Figure 5. 200 m temperature curves of different models ((a) uplift structure model; (b) fault structure model).

Table 2

Table 2. Calculation the parameter values of contribution rate and sensitivity coefficient in different models.

Table 3

Table 3. Contribution rates and sensitivity coefficients of different structures to temperature fields at different depths under the heat-flow coupling mode.

By comparing the contribution rate of middle-deep structure to temperature field at different depths, it can be found that the contribution rate coefficient increases with the increase of temperature field depth, and the maximum contribution rate coefficient of middle-deep structure to temperature field at 300 m depth is about 3 times that of temperature field at 100 m depth (Figure 6; Table 3). Comparing the contribution rate coefficients of different structures to the temperature field at the same depth, it can be found that the contribution rate of the fault structure to the temperature field is much larger than that of the uplift structure, generally remaining at around 5 to 8 times. Comparing the sensitivity coefficients of the temperature field at different depths to the middle and deep structures, it can be found that the sensitivity coefficient of the structure gradually increases with the increase of the depth of the temperature field. Comparing the sensitivity of different structures to temperature field, it can be found that the sensitivity of temperature field to fault structure is higher than that of uplift structure. The sensitivity coefficient of the fault structure is 2.4 times, 1.8 times and 1.5 times of the uplift structure at the depth of 100 m, 200 m and 300 m, respectively, and the average is 1.9 times (Table 3). Whether it is the contribution rate of the middle-deep structure structural formations to the temperature fields at different depths, or the sensitivity of the temperature fields to the middle and deep structural formations, the fault structures are significantly more prominent than the uplift structures. This is mainly because the fluid conduction in the fault structures plays a significant role, transferring the heat from the middle and deep layers to the shallow layers.

Figure 6

Figure 6. Contribution coefficient of different structures to shallow temperature field (a) 100 m deep temperature field; (b) 200 m deep temperature field; (c) 300 m deep temperature field.

Comparative analysis of different coupling fields

In order to better analyze the influence of middle-deep structure on shallow geothermal field, the thermo-solid coupling model is selected for simulation, and the simulation results are compared with the simulation results of the heat flow coupling model (Figure 7).

Figure 7

Figure 7. Thermo-solid coupling simulation results (a) Temperature field of the reference model; (b) Temperature field of uplift structure (buried depth of 2.1 km); (c) Temperature field of uplift structure (buried depth of 1.4 km); (d) Temperature field of fault structure temperature field (cutting depth 2.45 km); (e) Temperature field of fault structure temperature field (cutting depth 3.4 km); (f) Temperature field of fault structure temperature field (cutting depth 3.8 km).

Based on Equations 14, 15, the contribution rate and sensitivity coefficient were calculated. Contribution rate and sensitivity coefficient of different structures to temperature field at different depths under thermal-solid coupling mode are shown in Table 4. It can be seen that the average contribution rate to the temperature field increases with the increase of depth, whether it is a uplift structure or a fault structure. The average contribution rate of the uplift structure 1 to the temperature field is 0.35%, and the contribution rate of the uplift structure 2 to the temperature field is 4.87%.It can be seen that the shallower the buried depth of the uplift is, the more intense its influence on the shallow temperature field will be. The average contribution rate of the fault structure 1 to the temperature field is 22%, that of the fault structure 2 is 38.67%, and that of the fault structure 3 is 47%. It can be seen that the influence of fault structures on the temperature field is much stronger than that of uplift structures, and as the depth of the fault increases, the influence also gradually increases. Compared with fault structures, the temperature field is more sensitive to uplift structures.

Table 4

Table 4. Contribution rate and sensitivity coefficient of different structures to temperature field at different depths under thermal-solid coupling mode.

The average contribution rate of different structural models to the temperature field are compared (Figure 8). The average contribution rate of the uplift structure 1 to temperature field of thermal-solid coupling field is greater than that of heat flow coupling field (Figure 8a). However, the average contribution rate of uplift structure 2 to the temperature field of the thermal-solid coupling field is greater than that of the heat flow coupling field. This is because when the buried depth of the uplift structure is shallow, the heat transfer of rock is dominant, while when the buried depth of the uplift structure is deep, the heat transfer of water flow is dominant. By comparing the contribution rate of the fault structure, it can be found that the contribution rate of the fault structure to the shallow temperature field in the heat-flow coupling simulation mode is about 13 times that of the thermal-solid coupling (Figure 8b).

Figure 8

Figure 8. The average contribution rate coefficient of different structural models to the temperature field (a) Comparison of the contribution rate of heat flow coupling and thermal-solid coupling to the geothermal field in the uplift structure; (b) Comparison of the contribution rate of heat flow coupling and thermal-solid coupling to the geothermal field in the fault structure.

By comparing the sensitivity of the formation temperature field to the middle and deep structures under the two simulation modes of thermal-solid coupling and heat flow coupling, it can be found that sensitivity coefficient in the thermal-solid coupling model are significantly smaller than those in the heat flow coupling model. Additionally, it is found that the sensitivity coefficient of the temperature field to the fault structure in the thermal-solid coupling model is not only smaller than that of the heat-flow coupling simulation, but also much smaller than that of the uplift structure (Figure 9; Table 4). The seepage field in the stratum plays a great role in the heat transfer at the middle-deep layer, especially in the presence of fault structure.

Figure 9

Figure 9. Maximum sensitivity coefficient of different structures to shallow temperature field.

Conclusion

By constructing different middle-deep structure (uplift and fault structure) models, and based on numerical simulation methods such as thermal-solid coupling and heat-flow coupling, the general law of the influence of different middle-deep structures on shallow geothermal energy is simulated and analyzed. The contribution rate of the middle-deep structure to the shallow geothermal field and the sensitivity of the shallow geothermal field at different depths to the middle-deep structure are studied and analyzed. The findings are as follows.

1. Through the study of heat-flow coupling simulation, it is found that the contribution rate coefficient increases with the increase of temperature field depth, and the maximum contribution rate coefficient of the middle-deep structure to 300 m deep temperature field is about 3 times that of 100 m deep temperature field. The contribution rate of fault structure to temperature field is much greater than that of the uplift structure, generally remaining at around 5 to 8 times; As the depth of the temperature field increases, the sensitivity coefficient of the structure gradually increases; the sensitivity of the temperature field to the fault structure is higher than that of the uplift structure, and the sensitivity of the fault structure to the temperature field is about 1.9 times that of the uplift structure.

2. The contribution rate of the fault structure to the shallow temperature field in the heat-flow coupling mode is about 13 times that of the thermal-solid coupling. The sensitivity of the temperature field in the thermal-solid coupling model to the middle-deep structure is significantly lower than that in the heat-flow coupling model. The shallower the buried depth of the uplift is, the more intense its influence on the shallow temperature field will be.

3. The influence of fault structures on the temperature field is much stronger than that of uplift structures. This is mainly due to the fact that the fluid conduction in the fault structure plays a greater role in transferring the heat from the middle-deep layer to the shallow layer. Thermal convection plays an important role than thermal conduction in the influence of middle-deep structure on shallow geothermal fields.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

FG: Writing – review and editing, Conceptualization. JG: Writing – original draft, Formal Analysis. LW: Data curation, Writing – original draft. HB: Methodology, Writing – review and editing. XT: Funding acquisition, Project administration, Writing – review and editing. JQ: Project administration, Writing – review and editing. DZ: Conceptualization, Writing – original draft. XY: Investigation, Writing – original draft. TL: Investigation, Writing – original draft. QL: Data curation, Writing – original draft.

Funding

The authors declare that financial support was received for the research and/or publication of this article. This work was jointly funded by the Key Laboratory of Shallow Geothermal Energy, Ministry of Natural Resources of the People’s Republic of China (No. KLSGE202502-01 and No. KLSGE202402-02), China Geological Survey Project (No. DD20221677) and Technology Transformation Project (No. HX2025-13).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ahmed, A. A., Assadi, M., Kalantar, A., Sliwa, T., and Sapińska-liwa, A. (2022). A critical review on the use of shallow geothermal energy systems for heating and cooling purposes. Energies 15 (12), 4281. doi:10.3390/en15124281

CrossRef Full Text | Google Scholar

Chen, J. L., Luo, W. X., Dou, B., Zhou, Y., and Ning, W. T. (2021). Numerical simulation of geothermal field in a three-dimensional multi-fractured geological model of zhuolu basin. Bulletion Geol. Sci. Technol. 40 (3), 22–33. doi:10.19509/j.cnki.dzkq.2021.0317

CrossRef Full Text | Google Scholar

COMSOL Inc. (2020). COMSOL 5.6 user's guide. Burlington, MA: COMSOL Inc.

Google Scholar

Dai, M. G., Sun, P. G., Lei, H. F., Xing, Q., and Bao, Z. D. (2023). The spatial structure characteristics of strata and main thermal reservoirs in Xiong 'an New Area and the potential of geothermal water resources. Geol. Sci. 58 (2), 26. doi:10.12017/dzkx.2023.026

CrossRef Full Text | Google Scholar

Darius, M., R, P., and Vogt, C. (2011). The geothermal project den haag: 3D numerical models for temperature prediction and reservoir simulation. Geothermics 40 (3), 199–210. doi:10.1016/j.geothermics.2011.07.001

CrossRef Full Text | Google Scholar

Feng, H. W., Xu, H., Feng, H. T., and Gao, Y. (2024). Numerical simulation and experimental study of medium and deep ground source heat pump system in a cold and arid region. Case Stud. Therm. Eng. doi:10.1016/j.csite.2024.105473

CrossRef Full Text | Google Scholar

Figueira, J. S., Garcia, G. A., Vieira, A., Michopoulos, A. K., Boon, D. P., Loveridge, F., et al. (2024). Shallow geothermal energy systems for district heating and cooling networks: review and technological progression through case studies. Renew. Energy, 236. doi:10.1016/j.renene.2024.121436

CrossRef Full Text | Google Scholar

Franco, A., and Vaccaro, M. (2014). Numerical simulation of geothermal reservoirs for the sustainable design of energy plants: a review. Renew. Sustain. Energy Rev. 30, 987–1002. doi:10.1016/j.rser.2013.11.041

CrossRef Full Text | Google Scholar

Gao, J., Li, Y. Y., Wang, G. L., Zhang, B. J., and Xing, Y. F. (2023). Geothermal reservoirs characteristics and resource assessment of jixian system in gaoyang geothermal field, xiongan new area. Acta Geol. Sin. 44 (1), 133–144. doi:10.3975/cagsb.2022.100902

CrossRef Full Text | Google Scholar

Gude, V. G. (2016). Geothermal source potential for water desalination - current status and future perspective. Renew. Sustain. Energy Rev. 57 (May), 1038–1065. doi:10.1016/j.rser.2015.12.186

CrossRef Full Text | Google Scholar

Hu, Q. Y., Gao, J., Ma, F., Zhao, Z. H., Liu, G. H., Wang, G. L., et al. (2020). Dynamic prediction of recoverable geothermal resources in Rongcheng uplift area of Xiong'an new area. J. Geol. 94 (7), 2013–2025. doi:10.19762/j.cnki.dizhixuebao.2020230

CrossRef Full Text | Google Scholar

Ioan, S., and Alexandru, D. (2024). Investigation by monitoring and numerical simulation of the performance of a ground-source heat pump system with and without regeneration in different configurations and operating modes. Appl. Therm. Eng. 250, 1359–4311. doi:10.1016/j.applthermaleng.2024.123507

CrossRef Full Text | Google Scholar

José, A. A., José, R. C., Claudia, S. H., Jorge, R., and Thomas, G. K. (2022). Updated conceptual and numerical model of the los humeros geothermal field. Geothermics 106, 102564. doi:10.1016/j.geothermics.2022.102564

CrossRef Full Text | Google Scholar

Lei, X. D., Hu, S. B., Yang, Q. H., Jiang, G. Z., Li, J., Li, C., et al. (2017). Characteristics and hydrothermal system and its origin of nankou-sunhe fault zone in Beijing. Chin. J. Geophys. 60 (5), 1838–1850. doi:10.6038/cjg20170519

CrossRef Full Text | Google Scholar

Lei, X., Zheng, X., Duan, C. Y., Ye, J. H., and Liu, K. (2019). Three-dimensional numerical simulation of geothermal field of buried pipe group coupled with heat and permeable groundwater. Energies 12 (19), 3698–16. doi:10.3390/en12193698

CrossRef Full Text | Google Scholar

Li, W. (2021). Geothermal field in Tongzhou District of Beijing and its influencing factors. Geol. Bull. China 40 (7), 1189–1194. doi:10.12097/gbc.dztb-40-7-1189

CrossRef Full Text | Google Scholar

Li, N. B., Li, X., Yang, J. W., and Yu, Y. (2017). How to implement' shallow geothermal energy + ' model. Urban Geol. 12 (03), 7–13. doi:10.3969/j.issn.1007-1903.2017.03.002

CrossRef Full Text | Google Scholar

Li, K. F., Zhu, C. Q., Ma, Z., Yang, J. S., Su, H., and Xing, S. D. (2024). Thermophysical characteristics of rocks in Hainan Island and its influence on the geothermal field. Acta Geol. Sin., 1–21. doi:10.19762/j.cnki.dizhixuebao.2023085

CrossRef Full Text | Google Scholar

Liu, J., Zhang, W., Xue, Y., Wang, H. M., Li, S. T., Zhang, Y., et al. (2025). Evaluation of acid fracturing in carbonatite geothermal reservoirs based on a coupled thermo-hydro-mechanical-chemical model considering discrete fracture networks. Geomechanics Energy Environ. 43 (000), 100704. doi:10.1016/j.gete.2025.100704

CrossRef Full Text | Google Scholar

Luan, Y. B., Zheng, G. S., and Wei, W. S. (2013). Review of the shallow geothermal energy resources development and utilizaiton. Geol. Explor. 49 (02), 379–383.

Google Scholar

Lund, J. W., and Toth, A. N. (2021). Direct utilization of geothermal energy 2020 worldwide review. Geothermics 90, 101915. doi:10.1016/j.geothermics.2020.101915

CrossRef Full Text | Google Scholar

Ma, F., Wang, G. L., Zhang, W., Zhu, X., Zhan, H. X., and Yue, G. F. (2020). Structure of geothermal reservoirs and resource potential in the Rongcheng geothermal field in Xiong’an New Area. Acta Geol. Sin. 94 (7), 1981–1990. doi:10.19762/j.cnki.dizhixuebao.2020217

CrossRef Full Text | Google Scholar

Ma, F., Gao, J., Wang, G. L., Liu, G. H., Yu, M. X., Zhao, Z. H., et al. (2023). Numerical simulation of production and reinjection in carbonate geothermal reservoir of Rongcheng geothermal field, Xiong'an new area. J. Jilin Univ. 53 (5), 1534–1548. doi:10.13278/j.cnki.jjuese.20220004

CrossRef Full Text | Google Scholar

Pang, Z. H., Kong, Y. L., Shao, H. B., and Kolditz, O. (2018). Progress and perspectives of geothermal energy studies in China: from shallow to deep systems. Environ. Earth Sci. 77, 580. doi:10.1007/s12665-018-7757-z

CrossRef Full Text | Google Scholar

Roka, R. B., Figueiredo, A. J. P. D., Vieira, A. M. C. P. V., and Cardoso, J. C. D. P. C. (2023). A systematic review on shallow geothermal energy system: a light into six major barriers. Soils Rocks. 46 (1), e2023007622. doi:10.28927/SR.2023.007622

CrossRef Full Text | Google Scholar

Shi, Y. Z., Wang, S. J., Xiao, H. P., She, G. L., Rao, S., Hu, S. B., et al. (2022). 3D GeoModeller-based simulation of the geothermal field in the northern Songliao Basin. Nat. Gas. Ind. 42 (4), 46–53. doi:10.3787/j.issn.1000-0976.2022.04.004

CrossRef Full Text | Google Scholar

Sun, C. S., Shi, C. Y., Zhu, Z. M., Lin, H. X., Li, Z. H., Du, F., et al. (2025).Overburden failure characteristics and fracture evolution rule under repeated mining with multiple key strata control. Sci. Rep. 15,28029. doi:10.1038/s41598-025-14068-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang (2019). Study on monitoring of shallow ground temperature field and influencing factors of regional geothermal field in Nanjing city. Nanjing: Nanjing University, Master's thesis.

Google Scholar

Teng, T., Chen, Y. L., Wang, Y. M., and Qiao, X. X. (2025). In situ nuclear magnetic resonance observation of pore fractures and permeability evolution in rock and coal under triaxial compression. J. Energy Eng. 151 (4), 04025036. doi:10.1061/JLEED9.EYENG-6054

CrossRef Full Text | Google Scholar

Torresan, F., Piccinini, L., Cacace, M., Pola, M., Zampieri, D., and Fabbri, P. (2022). Numerical modeling as a tool for evaluating the renewability of geothermal resources: the case study of the Euganean Geothermal System (NE Italy). Environ. Geochemistry Health 44 (7), 2135–2162. doi:10.1007/s10653-021-01028-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Urchueguia, J. F. (2016). Shallow geothermal and ambient heat technologies for renewable heating. Renew Heating and Cooling, 89–118. doi:10.1016/B978-1-78242-213-6.00005-9

CrossRef Full Text | Google Scholar

Walch, A., Li, X., Chambers, J., Mohajeri, N., Yilmaz, S., Patel, M., et al. (2022). Shallow geothermal energy potential for heating and cooling of buildings with regeneration under climate change scenarios. Energy 244, 123086. doi:10.1016/j.energy.2021.123086

CrossRef Full Text | Google Scholar

Wang, W. L., Wang, G. L., Zhu, X., and Liu, Z. M. (2017). Evaluation of the development conditions and potential of shallow geothermal energy utilization in China's provincial capital cities. China Geol. 44 (06), 1062–1073. doi:10.12029/gc20170602

CrossRef Full Text | Google Scholar

Wang, X. Y., Yao, H. C., Li, J., Wang, Y. F., and Zhu, Y. Z. (2019). Experimental and numerical investigation on heat transfer characteristics of ammonia thermosyhpons at shallow geothermal temperature. Int. J. Heat Mass Transf. 136, 1147–1159. doi:10.1016/j.ijheatmasstransfer.2019.03.080

CrossRef Full Text | Google Scholar

Wang, G. L., Gao, J., Zhang, B. J., Xing, Y. F., Zhang, W., and Ma, F. (2020). Study on the thermal storage characteristics of the Wumishan formation and huge capacity gerothermal well parameters in the Gaoyang low uplift area of Xiong’an New Area. Acta Geol. Sin. 94 (7), 1970–1980. doi:10.19762/j.cnki.dizhixuebao.2020235

CrossRef Full Text | Google Scholar

Wang, G. L., Ma, F., Zhang, W., Zhu, X., Yu, M. X., Zhang, H. X., et al. (2024). Study on the dominant heat transfer mechanism of buried hills in North China: taking Xiong 'an New Area as an example. Earth Sci. Front. 31 (6), 052–066. doi:10.13745/j.esf.sf.2024.7.10

CrossRef Full Text | Google Scholar

Wang, W. L., Duan, Y. J., Zhang, W., Zhu, X., Ma, F., and Wang, G. L. (2024). Control factors and guidelines for urban-scale shallow geothermal energy development based on control units—The case of the starting area of Xiong’an New Area. Earth Sci. Front. 31 (6), 158–172. doi:10.13745/j.esf.sf.2024.7.16

CrossRef Full Text | Google Scholar

Wang, Z. T., Gu, C. C., Lu, H. F., Wang, Y. B., Hu, S. B., Jiang, G. Z., et al. (2025). Influence of thrust nappe structure on the geothermal field in Huainan-Huaibei area. Earth Sci. Front. doi:10.13745/j.esf.sf.2025.3.44

CrossRef Full Text | Google Scholar

Wang, L., Zhu, L., Cao, Z., Liu, J., Xue, Y., Wang, P., et al. (2025). Thermo-mechanical degradation and fracture evolution in low-permeability coal subjected to cyclic heating–cryogenic cooling. Phys. Fluids 37 (8), 086617. doi:10.1063/5.0282266

CrossRef Full Text | Google Scholar

Wei, W. S., and Li, X. (2021). Analysis on the development strategy of shallow geothermal energy in China in the next 10-15 years. Urban Geol. 16 (01), 1–8. doi:10.3969/j.issn.1007-1903.2021.01.001

CrossRef Full Text | Google Scholar

Wei, W. S., Zheng, G. S., Luan, Y. B., and Yang, J. W. (2012). Temperature characteristics of the room temperature zone and genetic mechanism of shallow geothermal energy. Urban Geol. 7 (02), 1–5. doi:10.3969/j.issn.1007-1903.2012.02.001

CrossRef Full Text | Google Scholar

Wei, W. S., Li, N. B., Zheng, G. S., Luan, Y. B., Yang, J. W., Li, X., et al. (2020). Study on genetic mechanism and control conditions of shallow geothermal energy in China. Urban Geol. 15 (01), 1–8.

Google Scholar

Wu, Z. W., and Song, H. Z. (2010). Numerical simulation of thermal convection in the shallow geothermal field. Rock Soil Mech. 231 (4), 6. doi:10.3969/j.issn.1000-7598.2010.04.051

CrossRef Full Text | Google Scholar

Xu, Y. S., Wang, X. W., Shen, S. L., and Zhou, A. (2020). Distribution characteristics and utilization of shallow geothermal energy in China. Energy Build. 229 (2), 110479. doi:10.1016/j.enbuild.2020.110479

CrossRef Full Text | Google Scholar

Yuan, X. Q., Xiao, Y. M., Lian, J. Y., Li, G., and Ma, M. R. (2024). Establishment of a geothermal basic database within 1000 m in Chinese mainland. Heat. Vent. Air Cond. 54 (12), 46–50. doi:10.19991/j.hvac1971.2024.12.08

CrossRef Full Text | Google Scholar

Yue, G. F., Zhao, Z. W., Wang, W. L., Zhao, M., and Lin, W. J. (2017). Numerical simulation study on soil temperature field influenced by GSHP:a case study in Beijing. Chin. Sci. Technol. Pap. 12 (9), 1049–1053. doi:10.3969/j.issn.2095-2783.2017.09.014

CrossRef Full Text | Google Scholar

Zeng, Y., He, B., Tang, L., Wu, N., Song, J., and Zhao, Z. (2020). Numerical simulation of temperature field and pressure field of the fracture system at Zhangzhou geothermal field. Environ. Earth Sci. 79 (11), 262. doi:10.1007/s12665-020-09018-y

CrossRef Full Text | Google Scholar

Zhang, H. Z. (2022). Utilization of shallow layer geothermal resources should be the preferred option for carbon reduction in Shanghai. Shanghai Energy Conservation (01), 5–8. doi:10.13770/j.cnki.issn2095-705x.2022.01.002

CrossRef Full Text | Google Scholar

Zhu, X., Gao, Z., Chen, T., Wang, W., Lu, C., and Zhang, Q. (2022). Study on the thermophysical properties and influencing factors of regional surface shallow rock and soil in China. Front. Earth Sci. 10, 864548. doi:10.3389/feart.2022.864548

CrossRef Full Text | Google Scholar

Glossary

ρ Density (kg/m³)

c_p Specific heat capacity (J/(kg·K))

T Temperature (K)

t Time (s)

u The fluid velocity (m/s)

${\mathit{\lambda }}_{\mathit{e}\mathit{f}\mathit{f}}$ Effective mean conductivity (W/(m·K))

d_f Width of the fault zone (m)

Q Fluid mass (kg/(m³·s))

$\varnothing$ Porosity (%)

${\mathit{\rho }}_{\mathit{w}}$ Water density (kg/m³)

k Permeability (m²)

${\mathit{\mu }}_{\mathit{w}}$ The dynamic viscosity of water (Pa·s)

p Fluid pressure (Pa)

g The acceleration of gravity (m/s²)

${\varnothing }_{\mathbf{f}}$ Porosity of the fault zone (%)

${\mathit{k}}_{\mathit{f}}$ The fault permeability (m²)

${\mathit{f}}_{\mathit{f}}$ The fracture roughness coefficient (Dimensionless)

K Bulk Modulus (Mpa)

${\mathit{\alpha }}_{\mathit{T}}$ Thermal expansion coefficient (K⁻¹)

${\mathrm{\epsilon }}_{\mathrm{v}}$ Volumetric Strain (Dimensionless)

$\mathbf{D}$ The fourth-order elastic matrix tensor (Pa)

$\mathit{\lambda }$ Lamé constants (Pa)

$\mathit{G}$ Lamé constants (Pa)

$\mathit{u}$ The displacement vector (m)

${\mathit{\rho }}_{\mathit{a}\mathit{v}}$ The average density (Kg/m³)

${\mathit{T}}_{\mathit{r}\mathit{e}\mathit{f}}$ The reference temperature (K)

M The contribution rate parameter (%)

E The sensitivity coefficient (Dimensionless)

T₀ The temperature of the initial temperature field (°C)

Δt The shallow temperature change after adding or changing the uplift or fault characteristics (°C)

H₀ The parameter value (depth) of the uplift or fault in the initial state (Km)

Δh The change of the parameter (Km)