In geothermal energy development, thermal conductivity is a critical physical parameter that determines the characteristics of subsurface rocks and their geothermal reservoirs. Thermal conductivity measurements provide fundamental data for evaluating the thermal transport properties of geothermal resources. Under high-temperature conditions, the thermal conductivity of rocks typically undergoes changes due to the movement of heat carriers, variations in mineral composition, and structural alterations. Additionally, rocks are often subjected to surrounding pressure and stress, which can modify their pore structure and, consequently, their thermal conductivity. This is particularly relevant during geothermal extraction, where the stress state of the rocks may undergo significant changes (Lin et al., 2025).

We have developed a high-temperature and high-pressure thermal conductivity instrument for hot dry rock. The instrument employs the steady-state heat flow meter method to measure the thermal conductivity of the material. The designed technical specifications for confining pressure and temperature are 150 MPa and 200°C, respectively. Granite samples from the hot dry rock reservoir in the Gonghe Basin, Qinghai, China, were tested using this instrument under both room temperature and high-temperature, high-pressure conditions. The Gonghe Basin hot dry rock (HDR) reservoir targets depths of 3,000–4,000 m, where lithostatic pressure ranges 70–100 MPa. The basin-specific thermal gradient is 45°C/km (Zhang et al., 2020; Song et al., 2021). Experimental conditions including temperature and pressure (80 MPa, 150°C) referred from borehole log information.

The measured thermal conductivities were 1.99 W/(m·°C) at room temperature and 1.69 W/(m·°C) under high-temperature and high-pressure conditions (Li S.-D. et al., 2024). These results indicate that the thermal conductivity of granite decreases by approximately 15% under high-temperature and high-pressure conditions compared to room temperature. The thermal conductivity of samples was utilized in the subsequent sections.

The heat extraction system for deep geothermal energy using cluster multi-branch U-shaped wells (UMW) involves coupled heat conduction between reservoir, wellbore and working fluid, and the design of complex well configurations is a typical three-dimensional problem. Therefore, it is necessary to develop a three-dimensional numerical algorithm for thermal-hydraulic coupling in reservoir-wellbore-working fluid.

We utilized the simulator allowing for three-dimensional numerical simulations of thermal-hydraulic coupling in deep geothermal reservoirs (https://gitee.com/geomech/hydrate) and undergone extensive benchmarking and validation at multiple scales (Xu et al., 2021; Li et al., 2023; Zhang et al., 2023; 2024b; 2024c).

In the numerical simulation of deep geothermal systems, the fundamental governing equations include the heat exchange between the working fluid in the wellbore and the reservoir. The energy conservation equation for the reservoir is expressed as follows (Equation 1): ρ r c r ∂ T r ∂ t + ρ r c r u · ∇ T r = ∇ · k r ∇ T r + q , (1) where ρ _r is the density of the reservoir rock (kg/m³), c _r is the specific heat capacity of the reservoir rock (J/(kg·°C)), T _r is the temperature of the reservoir rock (°C), k _r is the thermal conductivity of the reservoir rock (W/(m·°C)), t is time (s), q is the heat exchange term between the reservoir and the working fluid in the wellbore (W/m³).

For the heat transfer fluid in the wellbore, the energy conservation equation is expressed as (Equation 2): ρ f c f ∂ T f ∂ t + ρ f c f u · ∇ T f = ∇ · k f ∇ T f + q , (2) where ρ _f is the density of the working fluid (kg/m³), c _f is the specific heat capacity of the reservoir rock (J/(kg·°C)), T _f is the temperature of the working fluid (°C), k _f is the thermal conductivity of the working fluid (W/(m·°C)), t is time (s), u is the velocity vector of the working fluid (m/s).

The heat extraction power P can also be calculated based on the change in energy within the reservoir over time. The specific calculation formula is as follows (Equation 3): P = Δ E Δ t , (3) where ΔE is the change in energy within the reservoir (J), Δt is the time interval over which the energy change is measured (s).

The change in energy ΔE within the reservoir can be expressed as (Equation 4): Δ E = ∫ V ρ r c r Δ T r   d V , (4) where ρ _r is the density of the reservoir rock (kg/m³), c _r is the specific heat capacity of the reservoir rock (J/(kg·°C)), ΔT _r is the temperature of the reservoir rock (°C), V is the volume of the reservoir (m³).

In this model, heat conduction between the reservoir and the wellbore is prioritized, followed by heat conduction and convection processes of the working fluid within the wellbore. Heat transfer between the reservoir and the wellbore is exclusively governed by conduction, with no convective processes assumed within the reservoir. The temperature distributions of the reservoir and wellbore serve as mutual boundary conditions during the simulation, with data exchange occurring at regular intervals to ensure accurate evaluation of the heat exchange efficiency over long horizontal sections (Figure 2).

To validate the effectiveness and accuracy of the numerical method, we established a one-dimensional unsteady heat transfer model with a domain size of 100 m, discretized into 500 uniform grids. The left boundary was set as a fixed temperature boundary at 100°C, while the initial temperature of the model was 0°C. The thermal conductivity, density, and specific heat capacity of the model were configured based on experimental results from samples of the Gonghe Basin, with values set to 1.69 W/(m·°C), 2,640 kg/m³, and 754.4 J/(kg·°C), respectively. Temperature distribution results were extracted at six time points with intervals of 500 days. The numerical simulation results demonstrated excellent agreement with theoretical solutions, thereby validating the accuracy of the heat transfer simulation program to a certain extent (Figure 3a). The specific error analysis shows that the simulation results are close to analytic calculations. The mean absolute error (MAE) and root mean squared error (RMSE) are lower than 0.13°C and 0.3°C respectively (Figure 3b).

Traditional two-dimensional numerical models for geothermal development often simplify the wellbore direction through linear scaling, which fails to accurately reflect the actual heat exchange efficiency of long-distance horizontal wells. To address this limitation, we established a three-dimensional reservoir-working fluid thermal-hydraulic coupled model at the field scale, based on the geothermal distribution characteristics of the Gonghe Basin in Qinghai Province, China (Figure 4a). This model incorporates a U-shaped multi-lateral well system with six branch wells, enabling a more effective investigation of the influence of wellbore flow dynamics on horizontal well heat exchange efficiency and the calculation of effective power during the heat exchange process. Based on high-temperature and high-pressure thermal conductivity measurements of dry hot rock samples, the thermal conductivity of the reservoir in the numerical simulation is set to 1.69 W/(m·K), with a density of 2,640 kg/m³ and a specific heat capacity of 754.4 J/(kg·K). Additional parameters for the numerical model are provided in Table 1 (Song et al., 2021; Zhang et al., 2021).

The model employs a three-dimensional reservoir with dimensions of 3,000 m × 500 m × 200 m (Figure 4b). The model includes six horizontal branch wells arranged in an axisymmetric distribution, each with a length of approximately 2,500 m (Figure 4c). The reservoir model is discretized into a uniform grid with cell sizes of 15 m × 5 m × 4 m (Figure 4d).

The boundaries of the model are treated as adiabatic, and the initial reservoir temperature is set to 200°C, with the injected fluid temperature fixed at 50°C. To enhance computational efficiency and facilitate result visualization, only half of the actual geological model is simulated. The simulation period is set to 50 years.

Figure 5 presents the temporal evolution of the temperature field in a geothermal reservoir over 1 year, 5 years, and 50 years, highlighting the interaction between the injection well (blue rectangle) and the production well (red rectangle). At the 1-year mark, the injection well creates a localized cooling zone, as evidenced by the blue region surrounding it. The thermal front propagates outward along the flow path toward production well, forming distinct temperature gradients. This stage reflects the initial response of the reservoir to fluid injection and heat extraction. By the 5-year mark, the cooling front expands significantly along the flow channel, driven by advective heat transfer. While noticeable cooling occurs near the wells, regions farther from the flow path retain substantial thermal energy, showcasing the anisotropic nature of heat transfer within the reservoir. The combined effects of conduction and advection lead to progressive thermal changes over time.

Figure 6 provides a quantitative depiction of the reservoir thermal response over a 50-year operational period under a constant injection rate of 30 m³/h, tracking the three key indicators. The average reservoir temperature (blue line, circle markers, left y-axis), initiates at the original reservoir temperature of 200°C and exhibits a consistent, near-linear decrease throughout the simulation, reaching approximately 195.5°C by the 50-year mark. Concomitantly, the average temperature drop (orange line, square markers, right y-axis), which quantifies the deviation from the initial 200°C, shows a corresponding and steady increase from 0°C at the start to approximately 4.5°C after 50 years of operation. Complementing these temperature metrics, the proportion of grid cells experiencing a temperature decrease greater than 0.1°C (green line, diamond markers, far-right y-axis) demonstrates a progressive expansion of the thermal response zone, rising from 0% at the beginning of the simulation to approximately 45% by the end of the 50-year period. Collectively, these trends illustrate the continuous and progressive thermal drawdown and expansion of the thermally affected zone within the reservoir due to sustained heat extraction at the specified injection rate.

At 50 years, the reservoir undergoes pronounced thermal depletion along the primary flow pathway, with the cooling front reaching the production well. The extensive cooling zone near the injection well indicates significant heat extraction from the surrounding rock. The thermal breakthrough of injected cooler fluids at the production well reduces heat extraction efficiency, marking a critical point in the reservoir’s thermal lifespan. This long-term behavior underscores the balance between advective and conductive heat transfer mechanisms, as well as the impact of flow dynamics and reservoir heterogeneity. The analysis highlights the importance of optimizing injection and production strategies to delay thermal breakthrough, mitigate reservoir depletion, and ensure the sustainable utilization of geothermal energy resources over extended operational periods.

Figure 7 illustrates the spatial and temporal evolution of fluid temperature within six branch wells (a-f) over a 50-year operational period. The temperature distribution ranges from 50°C (blue) to 200°C (red), with yellow contour lines highlighting two critical temperature thresholds: ∼129°C, the normal operating temperature for the organic Rankine cycle (ORC), and ∼100°C, the minimum temperature required for ORC operation (Altun and Kilic, 2020). The analysis demonstrates distinct variations in thermal performance among the six wells due to differences in flow dynamics, heat extraction rates, and reservoir heterogeneity.

Over time, all six wells exhibit a downward trend in fluid temperature, with the 129°C and 100°C contours progressively moving toward the wellhead. The results for well 1 and well 6, well 2 and well 5, as well as well 3 and well 4, show identical thermal behaviors due to axisymmetric of six branch wells. Wells 1, 2, and 3 experience faster cooling, as the 129°C contour retreats rapidly, indicating earlier thermal breakthrough of cooler fluids. The movement of the 129°C contour is crucial for assessing the long-term efficiency of ORC power generation. As 129°C contour retreats, the area of the reservoir capable of sustaining ORC at its normal operating temperature diminishes, reducing system efficiency. This behavior suggests higher flow rates or more efficient heat extraction in these wells, leading to faster reservoir thermal depletion. The 100°C contour also moves closer to the wellhead over time, particularly in wells 1 and 3, further limiting the usable thermal energy in the reservoir, which suggests that within the cluster of branched wells, the presence of numerous heat-extracted branch wells results in a reduction of the average heat exchange capacity per well for a given total heat exchange area. Consequently, the temperature of the wellbore in the internal branched wells decreases more rapidly compared to that in the external branched wells (well 1 and well 6). The accelerated cooling of internal wells is fundamentally attributed to the overlapping thermal influence zones between closely spaced internal wells. This spatial configuration creates shared thermal depletion regions where cumulative heat extraction exceeds the natural thermal recharge capacity of the localized reservoir volume. According to Fourier’s law of heat conduction, the thermal gradient between adjacent depletion zones becomes attenuated, reducing the effective heat flux toward internal wells. Concurrently, flow dynamics demonstrate that clustered wells create converging flow paths that preferentially channel cooler reinjected fluids toward the internal well network. This establishes a positive feedback loop: localized cooling increases fluid viscosity, enhancing convective heat transfer that further accelerates thermal drawdown. In contrast, external wells benefit from undisturbed peripheral reservoir sections where radial heat conduction from deeper formations maintains steeper thermal gradients.

We also calculated the outlet temperature of well 1-3 horizontal branch wells and the variation of heat extracted power (Figure 8). At the onset of operation, the outlet temperatures of wells 1-3 exceed 129°C, enabling a maximum total average power of 5.95 MW when the outlet temperature excess129°C. According to feedback from on-site construction experience, when the outlet temperature exceeds 129°C, the ORC power generation capacity may be approximately 10% of the heat extraction power (595 kW). This high initial power output is attributed to the optimal thermal gradient and efficient heat exchange conditions within the reservoir. The elevated outlet temperatures indicate a substantial thermal energy reserve, which is effectively harnessed during the early stages of operation. This phase is critical for establishing baseline performance metrics and understanding the initial energy extraction potential of the geothermal system.

Over time, the outlet temperatures of wells 1-3 exhibit a gradual decline. This decline is primarily driven by the continuous extraction of thermal energy from the reservoir, leading to a reduction in the available heat flux. The rate of temperature decrease is influenced by the temperature difference between the reservoir rocks and the working fluid. As the outlet temperature drops below 129°C but remains above 100°C, the total average power decreases to 4.68 MW. the ORC power generation capacity may be lower than ∼10% of the heat extraction power. We assume a conversion efficiency of 5% at this phase for ORC power generation (234 kW). This reduction reflects the diminishing thermal driving force available for energy conversion in the Organic Rankine Cycle (ORC) system.

When the outlet temperature falls below 100°C but remains above 80°C, the power output further declines to 3.91 MW. We assume a conversion efficiency of 1% at this phase for ORC power generation (39.1 kW). This phase underscores the importance of temperature management in maintaining efficient power generation. The decline in the outlet temperature below 100°C indicates a shift in the reservoir thermal dynamics, where the heat extraction process becomes less efficient, and alternative strategies, such as altering the number of branch wells in operation, may need to be considered. During the subsequent heat extraction process, the horizontal branch wells located near the central region have already exchanged the majority of the reservoir thermal energy, resulting in the reservoir temperature (200°C) closes to that of the injected working fluid (50°C). Consequently, it is advisable to design an economical operation plan for multi branch wells, such as deactivating the low-efficiency branch wells gradually, to minimize costs associated with inefficient energy extraction.

From the perspective of the average heat extraction power, the average heat extraction power of a single set of six branch wells over a 50-year operating cycle is ∼4.32 MW. As for the ORC generation, the average generating power of the first two decades (21.5 years) is ∼ 284.4 kW. Roughly five similar clusters of branch wells (six branch wells) would be needed to achieve 1.2 MW of generating power. The average generating power of 50 years would be relatively low, ∼ 144.6 kW. Roughly nine similar clusters of branch wells (six branch wells) would be needed to achieve 1.2 MW of generating power.

It can be observed that the results from previous two-dimensional horizontal well calculations, which did not account for flow processes within the wellbore, are less accurate compared to those obtained when incorporating both the flow of working fluid within the wellbore and heat exchange between the reservoir and the wellbore. This is because the heat exchange due to fluid flow within the wellbore is subject to periodic limitations and may also be influenced by convective effects within the wellbore, leading to a relatively lower heat extraction power. Based on the analysis of the wellbore outlet temperature, it is recommended that a staged utilization of geothermal energy be considered during the extraction process. During high-temperature periods, electricity generation should be prioritized, while during subsequent stable periods, the focus can shift to heating and other energy demands.

The analysis of wells 1-3 provides valuable insights into the temporal variation of outlet temperatures and their impact on heat extraction power. The observed trends emphasize the need for proactive reservoir management strategies, such as optimizing well spacing, working strategy for branched wells, or integrating hybrid energy systems, to sustain power output over extended periods. Future research should focus on developing predictive models to simulate long-term thermal behavior and exploring innovative technologies to improve the efficiency and sustainability of UMW method.