A numerical study of non-Darcy flow in EGS heat reservoirs during heat extraction

Abstract

Underground non-Darcy fluid flow has been observed and investigated for decades in the petroleum industry. It is deduced by analogy that the fluid flow in enhanced geothermal system (EGS) heat reservoirs may also be in the non-Darcy regime under some conditions. In this paper, a transient 3D model was presented, taking into consideration the non-Darcy fluid flow in EGS heat reservoirs, to simulate the EGS long-term heat extraction process. Then, the non-Darcy flow behavior in water- and supercritical CO₂ (SCCO₂)-based EGSs was simulated and discussed. It is found that non-Darcy effects decrease the mass flow rate of the fluid injected and reduce the heat extraction rate of EGS as a flow resistance in addition to the Darcy resistance which is imposed to the seepage flow in EGS heat reservoirs. Compared with the water-EGS, the SCCO₂-EGS are more prone to experiencing much stronger non-Darcy flow due to the much larger mobility of the SCCO₂. The non-Darcy flow in SCCO₂- EGSs may thus greatly reduce their heat extraction performance. Further, a criterion was analyzed and proposed to judge the onset of the non-Darcy flow in EGS heat reservoirs. The fluid flow rate and the initial thermal state of the reservoir were taken and the characteristic Forchheimer number of an EGS was calculated. If the calculated Forchheimer number is larger than 0.2, the fluid flow in EGS heat reservoirs experiences non-negligible non-Darcy flow characteristic.

Keywords

enhanced geothermal system / non-Darcy flow / heat extraction / Reynolds number / Forchheimer number

Introduction

Geothermal energy exploitation is increasingly moving toward unconventional hot dry rock (HDR) resources. The so-called enhanced or engineered geothermal system (EGS) aims at exploiting the heat from HDR, and has a tremendous potential for power generation [¹–³]. However, commercialization of EGS power plants has strict requirements on fluid circulation rates, the production temperature, and the sustainability of geothermal systems. Moreover, considerable research and development efforts are needed to resolve the related scientific and technical issues [³–⁶]. Proper descriptions of the seepage flow in EGS heat reservoirs are crucial to accurate modeling and design, and subsequently to the construction and operation of EGS power plants.

The fracture configuration and the network of fractures in EGS heat reservoirs are typically complicated, and thus, difficult to accurately describe in reservoir modeling. An effective way to numerically model the heat reservoir is to treat it as a porous medium. Previous research work mostly assumed Darcy fluid flow in EGS reservoirs [⁷–⁹], indicating a simple proportional relationship between the seepage velocity and the pressure gradient is assumed. However, the Darcy flow assumption is only valid when the pressure gradients and resultant flow velocities are sufficiently small, i.e. when the Reynolds number is much less than the critical Reynolds numbers that is 1 for the flow around a sphere. As the Reynolds number gradually increases, the nonlinearity appears in the velocity-pressure gradient relationship, which is the so-called non-Darcy flow. Forchheimer [⁴] deduced from experimental results, that an additional term in relation to the squared velocity is required to represent the non-Darcy flow effect, when the flow rates are high. In recent years, some studies have confirmed the existence of non-Darcy flow in EGS heat reservoirs, particularly within the region near the well bore. For example, Kohl et al. [¹⁰] have investigated non-Darcy flow behavior in fractured rock, based on the field tests of the EGS site at Soultz-sous-Forèts, France. Zhang and Xing [¹¹] have employed the Forchheimer equation and a nonlinear finite element model to study the non-Darcy effects in a geothermal reservoir. The results indicate that the Forchheimer flow in the near-well region led to very different fluid velocity and pressure scenarios in the well bores.

The physical properties of heat transfer fluids, such as density and dynamic viscosity, can significantly affect the fluid flow behavior. Water and supercritical CO₂ (SCCO₂) are the two candidate heat transfer fluids being used in EGSs. References [⁷,⁸,¹⁰] have compared water with SCCO₂ as the heat transfer fluid of EGS. Different heat extraction performance for water-based EGS and SCCO₂-based EGS can be ascribed mainly to the different fluid mobility (i.e., density divided by dynamic viscosity). Since SCCO₂ has a much larger mobility than water, it is expected that much stronger non-Darcy flow exists in SCCO₂-EGS reservoirs. Though some work [¹⁰,¹¹] have been conducted to study the non-Darcy effects on the seepage flow in water-EGS, there is generally lack of a thorough and comprehensive study; and for SCCO₂-EGS, to the best of the authors’ knowledge, no work has been carried out to deal with the possibly occurring non-Darcy flow in the heat reservoir during EGS heat extraction [⁷,¹²–¹⁵].

The present paper aims at numerically studying non-Darcy flow in EGS heat reservoirs. It focuses mainly on identifying the different non-Darcy flow behaviors in water- and SCCO₂- based EGS, hoping to propose a criterion that can be used to judge whether non-Darcy flow effects should be considered or not when modeling heat extraction from EGS.

Numerical model of heat extraction from EGS

The present paper aims to simulate the heat extraction process in an idealized EGS reservoir that consist of an injection well, one or more production wells, a heat reservoir, and surrounding impermeable hot rocks, as described in Refs. [⁹,¹⁶–¹⁸]. It is assumed that the heat reservoir is well fractured and can be treated as a porous medium. Surrounding the heat reservoir are the impermeable hot rocks, so that the fluid loss is not considered. Single phase flow is considered in the reservoir, i.e. the liquid phase in water-EGS and super critical phase in SCCO₂-EGS. Local thermal non-equilibrium is considered between the solid rock matrix and fluid flowing in the factures. The heat losses across the walls of the production/injection wells are ignored. The governing equations of the model are listed below [⁹,¹⁶–¹⁹].

Fluid mass continuity equation

(1){\displaystyle \frac{\partial (\epsilon {\rho }_{\mathrm{f}})}{\partial t}}+\nabla \cdot ({\rho }_{\mathrm{f}}\mathit{u})=\mathrm{0.}

Fluid momentum equation

(2){\displaystyle \frac{\partial ({\rho }_{\mathrm{f}}\mathit{u})}{\partial t}}+{\displaystyle \frac{{\rho }_{\mathrm{f}}\mathit{u}}{\epsilon }}\nabla \mathit{u}=-\nabla p+\nabla \cdot \mu \nabla \mathit{u}-{\displaystyle \frac{\mu }{K}}\mathit{u}+\epsilon {\rho }_{\mathrm{f}}\mathit{g}+{\mathit{S}}_{\mathrm{ND}}.

Energy equation for heat transport in the solid rock or rock matrix in the heat reservoir

(3){\displaystyle \frac{\partial [(\mathrm{1}-\epsilon ){\rho }_{\mathrm{s}}{c}_{p\mathrm{s}}{T}_{\mathrm{s}}]}{\partial t}}=\nabla \cdot ({\lambda }_{\mathrm{s}}^{\mathrm{eff}}\nabla T_{\mathrm{s}})-ha(T_{\mathrm{s}}-T_{\mathrm{f}}).

Energy equation for heat transport in the fluid

(4){\displaystyle \frac{\partial [\epsilon {\rho }_{\mathrm{f}}{c}_{p\mathrm{f}}{T}_{\mathrm{f}}]}{\partial t}}+\nabla \cdot (\mathit{u}{\rho }_{\mathrm{f}}{c}_{p\mathrm{f}}T)=\nabla \cdot ({\lambda }_{\mathrm{f}}^{\mathrm{eff}}\nabla T_{\mathrm{f}})+ha(T_{\mathrm{s}}-T_{\mathrm{f}})，

where u, p, T_f, and T_s are the primary variables to be solved, denoting the superficial fluid velocity vector, fluid pressure, fluid temperature, and rock temperature, respectively; r_f and r_s are the density of the fluid and solid rock, respectively; e and K denotes the porosity and permeability, respectively; c_pf and c_ps are the specific heat capacity of the fluid and solid rock, respectively; and l_f and l_s are the heat conductivity of the fluid and solid rock, respectively. The fluid properties in the above equations, including r_f, m, c_pf, and l_f, are considered to be dependent on the fluid pressure and temperature. The effective heat conductivity for the rock matrix, l_s^eff, and for the fluid, l_f^eff, are determined based on their volume fractions raised to the power of 1.5, i.e., l_s^eff = l_s(1 – e)^1.5 and l_f^eff = l_fe^1.5, respectively. The last terms, i.e., ±ha (T_s – T_f), in Eqs. (3) and (4), describes the heat exchange between the rock matrix and the fluid in the reservoir, where ha denotes the volumetric heat transfer coefficient. The so-called single-domain multi- sub-regions approach [⁹] is employed to treat the model geometry, namely: e = 0 and K = 0 in the solid rock region, e = 1 and K = ∞ in the injection and production wells, and e and K have finite values in the porous reservoir. This treatment of single-domain multi- sub-regions allows the unified mathematical description (Eqs. (1) and (2)) to fluid flow in different regions, simplifying the computational treatment of the interfaces between these regions.

The term S_ND in Eq. (2) is the so-called non-Darcy term, which is commonly described by the Forchheimer model [⁴,¹¹], that is

(5){\mathit{S}}_{\mathrm{ND}}=-\beta {\rho }_{\mathrm{f}}\left|\mathit{u}\right|\mathit{u},

where b is the non-Darcy coefficient. The non-Darcy term, which formulates the flow resistance caused by the liquid-solid interactions, has been generally adopted for describing the high-velocity flow behavior and validated experimentally and numerically by many researchers [¹⁰,¹¹,²⁰–²²]. The non-Darcy flow coefficient b, which is associated with the intrinsic characteristics of porous medium, is a key parameter to predict the flow behavior under non-Darcy flow regime. Generally, b is taken as a function of permeability and porosity or permeability only, which is usually determined from laboratory measurements and/or practical well-tests.

Case setup

EGS subsurface geometry and boundary conditions

The modeled EGS subsurface geometry consists of three sub-regions [⁹,¹⁶–¹⁹]: the injection and production wells (e = 1 and K = ∞), the geothermal reservoir (finite e and K), and the rock that enclosing the reservoir (e = 0 and K = 0). With this treatment of the EGS subsurface geometry, the interfaces between the wells, the geothermal reservoir, and the surrounding rock are interior interfaces and are automatically handled in the numerical modeling. Specially, the present paper considers a 5-spot well configuration EGS (one vertical injection well and four vertical production wells). The geometry and geometrical dimensions are schematically shown in Fig. 1(a). The simulated domain is a 2000 m × 6000 m × 2000 m volume, where the horizontal plane at z = 0 indicates the earth surface. The reservoir, a 500 m × 500 m × 500 m porous cubic region, is located at a depth of about 4000 m of the reservoir. The vertical injection well is positioned at the center of the horizontal plane and the four production wells are arranged close to the four corners. There is a 50 m distance between the production well center and the nearby heat reservoir boundary. The injection and production wells are all 0.2 m × 0.2 m square-shaped in the xy-plane. Only a quarter of the EGS geometry is simulated because of the geometrical symmetry, as shown in Fig. 1(b). The boundary surfaces, including the xz-plane with the minimum y value and the yz-plane with the minimum x value, are set as symmetry boundaries while the other external surfaces of the computational domain are specified as zero-flux boundaries for all primary variables, except for temperature, which has fixed values at these boundaries surfaces. The nonslip flow boundary condition is prescribed at the walls of the injection/production well boreholes. The fluid inlet and outlets have fixed pressure boundary conditions. The numerical mesh system is displayed in Fig. 1(b). The meshes are carefully designed to ensure sufficiently fine meshes in the injection and production wells and throughout the reservoir. The grid-independence tests performed indicate changes of less than 2% at the resultant pressure, velocity, and temperature distributions in the reservoir when the mesh system increases its element number from approximately 150000 to approximately 200000. Further increasing the numerical element number to 270000 generates negligible changes at the results. Therefore, the mesh system of 270000 elements is adopted for the numerical study.

Model parameters and cases simulated

The density, specific heat capacity, and thermal conductivity of the rock are set as r_s = 2650 kg/m³, C_ps = 1000 J/(kg^-1·K^-1), k_s = 2.4 W/(m ^-1·K^-1), respectively [⁹]. The geothermal gradient is assumed to be 4 K/100 m, while the mean annual surface temperature is fixed at 300 K [⁹]. Initially, the temperature of the subsurface working fluid is assumed to be locally the same as the rock temperature. The volumetric heat transfer coefficient ha is set as 1.0 W/(m ^-3·K^-1) [¹⁹]. During heat extraction, the cold fluid, with a temperature of 343.15 K, is injected into the geothermal reservoir at a prescribed pressure difference of 10 MPa between the injection and the production wells. The pressure at the center of the geothermal reservoir is set at 40 MPa.

Reservoir permeability and porosity may be the two most important parameters, as they dominate the EGS heat extraction process. They dictate the flow distribution and the flow resistance (i.e. the needed external pump work) and thereby, they directly affect the EGS performance, including the heat extraction performance, geothermal reservoir lifetime, and economic performance. It is assumed that the EGS reservoir has been homogeneously fractured, being of a constant and uniform porosity and permeability of 0.01 and 1.0 × 10^-14 m² [¹⁹], respectively. The correlations of the non-Darcy coefficient b are taken from Refs. [²⁰–²²]. Friedel and Voigt [²⁰] have summarized the experimental data from different sources and obtained a correlation for b. Janicek and Katz [²¹] have established a correlation for b as a function of permeability and porosity for natural porous media. Pascal et al. [²²] have described a method for predicting both the non-Darcy flow coefficient, b, and the vertical fracture length, based on a single-point, variable flow drawdown tests of shallow, low permeability gas reservoirs.

The b correlation is believed to be associated with the fracture properties in reservoirs. To compare the Darcy and non-Darcy flow behaviors and further investigate the criterion for judging the onset of non-Darcy flow, the above-mentioned three correlations of b have been considered in the numerical simulations. With K = 10.0 mD and e = 0.01, the b values are calculated in terms of the three correlations of b, respectively. The b correlations together with the calculated b values are listed in Table 1. Note that in Table 1, b has the unit of cm^-1, and K in mD. The b values calculated by different correlations show large differences.

Two groups of cases have been considered, as listed in Table 2, one group for the study of non-Darcy effects on EGS heat extraction process, and the other mainly for the derivation of a criterion that can be used to judge whether non-Darcy flow effects shall be considered or not when modeling EGS heat extraction processes.

The model equations, Eqs. (1–4), along with the above-described initial and boundary conditions, are solved with the commercial CFD software Fluent 6.3®. The equations of the mass continuity and the momentum conservation (Eqs. (1) and (2)) are solved by Fluent® itself and the two energy equations (Eqs. (3) and (4)) are solved by customizing the flexible, user-defined functions (UDFs). The SIMPLEC (semi-implicit method for pressure linked equation consistency) is used to address pressure-velocity coupling. The first-order upwind differencing scheme is generally used for the discretization of the spatial-derivative terms and a fully implicit scheme is used for the discretization of the transient terms. A variable time step is considered during the simulation. During the early period of EGS operations, a small time step of 1.0 s is implemented to ensure that the thermal front will not penetrate a cell during one single time step. After the thermal front passes the near-well region, i.e. after about two years, the time steps are progressively increased to 1.0 × 10⁷ s to accelerate computation.

Results and discussion

EGS heat extraction performance with or without non-Darcy effects

During the heat mining processes in EGS, the heat transfer fluid is driven by the pressure gradient in the heat reservoir, extracting the thermal energy from hot rocks. Figure 2 compares the simulated fluid pressure distribution in the reservoir at three representative time instants of 1, 5 and, 10 years into the EGS operation for Cases 1–4. The pressure values are relative to the reference pressure, 40 MPa, which is the assumed absolute pressure value at the center of the reservoir. In both water- and SCCO₂-based EGS, the fluid pressure remains high only in a very limited region close to the injection well borehole, indicating that most of the flow resistance originates from this region. Moreover, the pressure gradient is mainly in the horizontal direction, meaning that the fluid flow in the reservoir will mostly take place in the horizontal direction. The non-Darcy effect can bring more flow resistance and thus causes a larger pressure drop within a narrower region around the injection well. It can be seen that there is little difference at the pressure distributions of Cases 1 and 2, indicating slight non-Darcy effects on water-EGS. However, the locations with a high fluid pressure drop for Case 4 are constrained within a narrower region around the injection well compared with Case 3, indicating that the non-Darcy flow may cause more flow resistance in SCCO₂-EGS.

Figure 3 shows the seepage speed of Cases 1–4 at three representative time instants of 1, 5, and 10 years into the EGS operation. In both water and SCCO₂ based EGS, the seepage speed is higher in the vicinity region of the injection well and the production wells, where the wells have been approximated as a point source or sink. In the central region between the injection and production wells, the fluid speeds are relatively smaller. This low-velocity region first appears in the upper part of the reservoir and gradually broadens toward the bottom of the reservoir with the progress of EGS heat extraction. Two processes cause this flow behavior. One is that the seepage flow is faster in the deeper region of the reservoir, due to the larger fluid mobility, caused by higher fluid temperatures at greater depths, due to the imposed geothermal gradient. The other process is that the fluid seepage velocity reduces gradually along with the proceeding of heat extraction. For the water-EGS, very slight differences in seepage speed can be observed between Cases 1 and 2, indicating that the non-Darcy effect is insignificant. The seepage speed in SCCO₂-EGS is much higher than that in water-EGS, resulting in a more significant non-Darcy effect in SCCO₂-EGS (Fig. 3(b)).

Figure 4 depicts the temporal development of fluid temperature distribution in the heat reservoir of Cases 1–4. The injected cold fluid absorbs the heat stored in the rock matrix. A low-temperature region of the fluid is observed to appear around the injection well at first, which then gradually expands toward the production well as heat is extracted. Little difference in the fluid temperature field can be observed between Cases 1 and 2 (Fig. 4(a)). However, the propagation speed of the cold front, which denotes the edge of the low fluid temperature region, is reduced in SCCO₂-EGS, when the non-Darcy effect is considered. This means that the non-Darcy effect lowers the heat extraction rate in SCCO₂-EGS.

Figure 5 depicts the mass flow rate curves calculated for Cases 1–4. The non-Darcy effect has little impact on the mass flow rate in water-EGS but has a strong impact on that in SCCO₂-EGS. At the assumption of prevailing Darcy flow, the mass flow rate of SCCO₂-EGS is about 4 times that of water-EGS, at the same given identical pressure differences between injection and production wells. However, when the fluid flow is non-Darcian, the mass flow rate of the SCCO₂-EGS is only about three times that of the water-EGS. This means that the advantage of SCCO₂-EGS, i.e., significantly higher mass flow rates than in water-EGS, would be evidently weakened by the non-Darcy flow effect.

Figure 6 illustrates the production temperature curves calculated for Cases 1–4. The displayed EGS production temperature is the volumetrically averaged fluid temperature at the outlet of production wells. At the beginning of fluid production, the production temperature is about 460 K for all water- and SCCO₂- based cases, which corresponds to the initial average rock temperature in the reservoir. Thereafter, the production temperature decreases with time. The duration for the production temperature to remain at about 460 K, is longer for water-EGS cases, compared with SCCO₂-EGS cases. Moreover, this duration is longer if non-Darcy flow is assumed, compared with the cases in which only Darcy flow is considered. It can be observed that the non-Darcy effects have a stronger impact on the SCCO₂-EGS production temperature, compared with the water-EGS. By defining the EGS abandonment temperature as the EGS production temperature equaling 10 K below its maximum value, the EGS lifetime can be determined from Fig. 6 for the four cases. For the water-EGS, the lifetime is about 18.3 years under the Darcy flow assumption, while it is about 18.5 years when the non-Darcy flow is assumed. For the SCCO₂-EGS, the lifetime is about 11.8 years under the Darcy flow assumption, and about 13.1 years for the non-Darcy case. The longer lifetime of the non-Darcy case is mainly caused by the reduced fluid mass flow rate, as exhibited in Fig. 5.

Further, the real-time heat extraction rate was calculated, defined as the fluid enthalpy at the outlets of the production wells minus that at the inlet of the injection well, i.e., (Q_out(t)H_out(t)–Q_inj(t)H_inj). Figure 7 shows the non-Darcy effects on the real-time heat extraction rate of EGS. The temporal evolution of the heat extraction rate shows three distinct stages for all four cases. During the first stage, about 1 year into the EGS operation, the heat extraction rate suffers from a sharp reduction, due to the rapid decrease in fluid mass flow rate (see in Fig. 5). The second stage, with differing duration for each case, is characterized by a quasi-stable heat extraction rate. During the last stage, the heat extraction rate decreases at an accelerated rate, indicating that the EGS is close to the end of its lifetime. It is seen in Fig. 7 that, for water-EGS, the non-Darcy effect is insignificant whereas for SCCO₂-EGS the non-Darcy flow reduces the heat extraction rate, and, to some extent, also increases the duration of the stage, wherea quasi-stable heat extraction rate can be observed.

Criterion for judging the onset of non-Darcy flow

There are two main dimensionless numbers that may be used to indicate the onset of non-Darcy flow, namely the Reynolds number, Re, and the Forchheimer number, Fo. The Reynolds number for fluid flow in porous media is commonly defined similar to that for fluid flow in pipes. However, it is difficult to determine the characteristic length due to the often complicated structure of porous heat reservoirs. To facilitate its use, Котяхов [²³] has modified the definition of the Reynolds number for fluid flow in porous media as

(6)Re=\mathrm{4}\sqrt{\mathrm{2}}{\displaystyle \frac{{\rho }_{\mathrm{f}}\left|\mathit{u}\right|\sqrt{K}}{{\mu }_{\mathrm{f}}{\epsilon }^{\mathrm{3}/\mathrm{2}}}}.

It has been said that the seepage flow obeys Darcy’s law when Re is less than a critical value, Re_c = 0.2–0.3 [²³,²⁴].

The Forchheimer number is defined as [¹¹,²⁵]

(7)Fo={\displaystyle \frac{{\rho }_{\mathrm{f}}K\beta \left|\mathit{u}\right|}{{\mu }_{\mathrm{f}}}}.

The Forchheimer number was also applied to identify whether the non-Darcy behavior must be considered or not [¹¹,²⁵]. The critical value of the Forchheimer number is in the range of 0.005–0.2 [¹¹,²⁵]. For values of the Forchheimer number below the critical value, non-Darcy effects should be negligible. As shown by Eq. (7), the Fo depends on three factors: the mobility of the heat transfer fluid, r_f/m_f; the flow resistance properties, K and b, of the reservoir; and the seepage velocity, |u|. It is worth noticing that both the Reynolds number and the Forchheimer number are directly proportional to the mobility of the fluid, if the porosity and the permeability of the reservoir as well as the seepage velocity in the reservoir are constant.

Intuitively, the non-Darcy effect will impose an additional resistance to the fluid flow and thus reduce the flow velocity. To quantify the non-Darcy effect, a new parameter, g, namely, the relative velocity reduction, was defined as

(8)\mathit{\gamma }=\left|{\displaystyle \frac{U-U_{\mathrm{Darcy}}}{U_{\mathrm{Darcy}}}}\right|\times \mathrm{100}\%,

where U and U_Darcy denote the seepage velocity magnitude of the non-Darcy case and the corresponding Darcy case, respectively.

For the six cases (i.e., Cases 2, 4, 5–8, listed in Table 2) with non-Darcy flow considered, the g values was calculated. Defining a monitoring line (i.e., Line AC) in the reservoir, the g profile was plotted in Fig. 8 along this line for all the six cases at 5 years into the EGS operations. The position of Line AC is indicated in Fig. 2(a). It is located at the mid-xy plane of the heat reservoirs. Point A is positioned at the wall of the injection well borehole and Point C at the wall of the production well borehole. There is one more monitoring point, Point B, which is approximately at the mid-point of Line AC. The g values are larger in the near-regions of the injection/production wells. The position of the largest g value is always in the near-region of the production well. Except for Cases 2 and 5, where generally small g values exist in the reservoirs, the other four cases all yield relatively large g values (>10%), indicating more significant non-Darcy effects. The water-EGS cases do not show significant non-Darcy effects unless the non-Darcy flow coefficient, b, is sufficiently large. The SCCO₂-EGS cases all show significant non-Darcy effects.

To shed light on the correlation of non-Darcy effects and Re (or Fo), the calculated g, Re, and Fo values at the three monitoring points for the six cases (i.e., Cases 2, 4, 5–8) are tabulated in Table 3. It is indicated that increasing the non-Darcy coefficient, b, leads to more significant non-Darcy effects. At Point C, which is located at the wall of the production well borehole, g takes on the largest value for each case. For water-EGS, g is small, being 1.1% and 2.6% at Point C, if b takes the values of 1.3 × 10⁸ cm^-1 and 3.2 × 10⁸ cm^-1, respectively. However, when a 3.2 × 10⁹ cm^-1 non-Darcy coefficient is given, the maximum g value is observed to be 14.9% at Point C, and even at Point B, g is 10.1%. The g value is large for all the three SCCO₂-EGS cases considered, indicating significant non-Darcy effects.

Table 3 shows that the calculated Re values for water-EGS cases are small (<0.2), whereas for SCCO₂-EGS cases, they are large with a maximum of 2.25. There is no definite relation between the observed non-Darcy behavior and the Re values. Therefore, the Reynolds number is not a suitable parameter for judging the onset of non-Darcy flow in EGS heat reservoirs, which is an expected result, as the definition of the Reynolds number does not contain any information about the non-Darcy coefficient b.

Investigating the Fo values shown in Table 3, for Cases 4, and 6–8, which exhibit significant non-Darcy effects, the calculated Fo values are all large at Points A and C, and the larger the Fo, the more significant the non-Darcy effect is. For the other two cases (Cases 2 and 5) the calculated Fo values are small and no non-Darcy effects are observed. Therefore, the Forchheimer number has a better correlation with g and can thus be used as a criterion for judging the onset of the non-Darcy flow.

From Table 3, it is also observed that the local g does not positively correlate to the local Fo value. A larger local Fo value does not mean a larger g at this location, for example, Points A and C of Cases 7 and 8; a very small local Fo value does not mean the local g is negligible, for example, Point B Case 7 or Case 8. A better definition of Fo that can characterize the overall flow regime is critically needed.

The Forchheimer number evolves during EGS heat extraction operation and has a strong location-dependency. To calculate Fo, the local fluid mobility and velocity should be obtained. In the near-region of the injection well, the fluid experiences larger temperature changes (Fig. 4) and the fluid mobility may change significantly during EGS operations. Compared with Point A, determining the Fo at point C may be more robust to determine the onset of non-Darcy flow. Taking this into consideration, to facilitate practical uses, it is proposed in this paper to calculate the Forchheimer number at the wall of the production well in terms of the average fluid velocity and the fluid mobility at the initial average temperature in the reservoir. This Forchheimer number is then named as the characteristic Forchheimer number in EGS, Fo_ch, and defined as

(9)F{o}_{\mathrm{ch}}={\displaystyle \frac{QK\beta }{A{\mu }_{\mathrm{ref}}}},

where Q denotes the mass flow rate; A the cross-sectional surface area of the production well borehole, located within the heat reservoir, and m_ref the reference fluid dynamic viscosity, which determines in terms of the initial average reservoir temperature. The calculated Fo_ch values for the six cases are also listed in Table 3. Zhang et al. [¹¹] and Saboorian-Jooybari and Pourafshary [²⁵] have suggested that the critical Forchheimer number is within 0.005–0.2. A critical Forchheimer number of 0.2 is suggested in this paper, above which non-Darcy effects need to be considered.

It is worth pointing out that the Forchheimer model, Eq. (5), is only one description, though commonly used, of the non-Darcy flow. The other description is the Izbash model [²⁶], which uses a high-order polynomial equation of fluid velocity to calculate the non-Darcy flow resistance. However, its theoretical foundation is not yet well established [²⁷]. Considering the macroscopic shearing effect between the fluid and the porous matrix, Brinkman [²⁸] has added the 2nd-order derivative of fluid velocity to the Darcy model, resulting in the Brinkman equation, which is suitable to describe the non-Darcy flow in highly porous media. Using different non-Darcy flow models may lead to differing descriptions of the non-Darcy flow behavior and regimes. However, because the Forchheimer model has some definite physical meaning and is more widely applicable, it is recommend in this paper as a proper model for describing non-Darcy flow in EGS heat reservoirs.

Conclusions

To study the non-Darcy flow effects in EGS heat reservoir during heat extraction, a previous 3D, transient model has been extended in the present paper, by including additionally a non-Darcy Forchheimer resistance term in the momentum conservation equation. Simulations were performed for water- and SCCO₂- based EGS. The simulated cases are of an idealized 5-spot well layout, with one vertical injection well located at the center of the reservoir and four vertical production wells surrounding the injection well. Besides, the non-Darcy flow behaviors of water- and SCCO₂-based EGS have been compared, focusing mainly on the conditions that could result in non-Darcy flow.

Compared with water-EGS, SCCO₂-EGS is generally more likely to show a non-Darcy flow in heat reservoirs, as the SCCO₂ fluid has much larger mobility values. The non-Darcy effect leads to a reduction of the average fluid velocity, decreasing the mass flow rate of fluid injected and diminishing the heat extraction rate of EGS. The heat extraction performance of SCCO₂-EGS significantly deteriorates when non-Darcy flow is present in the heat reservoir. For the water-EGS cases, only when the non-Darcy flow coefficient is relatively large, is a non-Darcy effect detectable.

To establish a criterion that can be used to judge whether non-Darcy effects should be considered when modeling EGS heat extraction processes, two candidate parameters have been analyzed: the Reynolds number and the Forchheimer number. The Forchheimer number is found to have a definite correlation with the observed non-Darcy flow effects. At the locations where the Forchheimer number is large, the non-Darcy effects are more pronounced. Since the Forchheimer number is location-dependent, and varies with time during EGS heat extraction, to facilitate its practical uses, a characteristic Forchheimer number has been further defined, which can be calculated from the injection fluid mass flow rate and fluid properties at the initial average temperature in the heat reservoir. The critical value of the characteristic Forchheimer number is suggested to be 0.2, above which non-Darcy effects need to be considered when numerically modeling EGS.

References

[1] Tester J W, Anderson B J, Batchelor A S, Blackwell D D, Di Pippo R, Drake E M, Garnish J, Livesay B, Moore M C, Nichols K, Petty S, Toksoz M N, Veatch R W. The Future of Geothermal Energy: Impact of Enhanced Geothermal Systems (EGS) on the United States in the 21st Century. DOE Contract DE-AC07-05ID14517 Final Report, Massachusetts Institute of Technology, 2006

[2] Yang Y, Yeh H. Modeling heat extraction from hot dry rock in a multi-well system. Applied Thermal Engineering, 2009, 29(8–9): 1676–1681

[3] Olasolo P, Juárez M C, Morales M P, D´Amico S, Liarte I A. Enhanced geothermal systems (EGS): a review. Renewable & Sustainable Energy Reviews, 2016, 56: 133–144

[4] Forchheimer P. Water moving through soil. Zeitschrift des Vereines Deutscher Ingenieuer, 1901, 45: 1782–1788

[5] Gelet R, Loret B, Khalili N. A thermo-hydro-mechanical coupled model in local thermal non-equilibrium for fractured HDR reservoir with double porosity. Journal of Geophysical Research, 2012, 117: B07205

[6] Held S, Genter A, Kohl T, Kölbel T, Sausse J, Schoenball M. Economic evaluation of geothermal reservoir performance through modeling the complexity of the operating EGS in Soultz-sous-Forêts. Geothermics, 2014, 51: 270–280

[7] Brown D W. A hot dry rock geothermal energy concept utilizing supercritical CO₂ instead of water. In: The 25th Workshop on Geothermal Reservoir Engineering, Stanford University, Stanford, CA, USA, 2000

[8] Pruess K. On production behavior of enhanced geothermal systems with CO₂ as working fluid. Applied Thermal Engineering, 2008, 49: 1446–1454

[9] Jiang F M, Luo L, Chen J L. A novel three-dimensional transient model for subsurface heat exchange in enhanced geothermal systems. International Communications in Heat and Mass Transfer, 2013, 41: 57–62

[10] Kohl T, Evans K F, Hopkirk R J, Jung R, Rybach L. Observation and simulation of non-Darcian flow transients in fractured rock. Water Resources Research, 1997, 33(3): 407–418

[11] Zhang J, Xing H L. Numerical modeling of non-Darcy flow in near-well region of a geothermal reservoir. Geothermics, 2012, 42: 78–86

[12] Haghshenas Fard M, Hooman K, Chua H T. Numerical simulation of a supercritical CO₂ geothermosiphon. International Communications in Heat and Mass Transfer, 2010, 37(10): 1447–1451

[13] Xu C S, Dowd P, Li Q. Carbon sequestration potential of the Habanero reservoir when carbon dioxide is used as the heat exchange fluid. Journal of Rock Mechanics and Geotechnical Engineering, 2016, 8: 50–59

[14] Borgia A B, Oldenburg C M, Zhang R, Pan L, Daley T M, Finsterle S, Ramakrishnan T S. Simulations of CO₂ injection into fractures and faults for improving their geophysical characterization at EGS sites. Geothermics, 2017, 69: 189–201

[15] Wang C L, Cheng W L, Nian Y L, Yang L, Han B B, Liu M H. Simulation of heat extraction from CO₂-based enhanced geothermal systems considering CO₂ sequestration. Energy, 2018, 142: 157–167

[16] Cao W J, Huang W B, Jiang F M. Numerical study on variable thermophysical properties of heat transfer fluid affecting EGS heat extraction. International Journal of Heat and Mass Transfer, 2016, 92: 1205–1217

[17] Chen J L, Jiang F M. A numerical study of EGS heat extraction process based on a thermal non-equilibrium model for heat transfer in subsurface porous heat reservoir. Heat and Mass Transfer, 2016, 52(2): 255–267

[18] Chen J L, Jiang F M. Designing multi-well layout for enhanced geothermal system to better exploit hot dry rock geothermal energy. Renewable Energy, 2015, 74: 37–48

[19] Jiang F M, Chen J L, Huang W B, Luo L. A three-dimensional transient model for EGS subsurface thermo-hydraulic process. Energy, 2014, 72: 300–310

[20] Friedel T, Voigt H D. Investigation of non-Darcy flow in tight-gas reservoirs with fractured wells. Journal of Petroleum Science Engineering, 2006, 54(3-4): 112–128

[21] Janicek J D, Katz D L. Applications of unsteady state gas flow calculations. In: Research Conference on Flow of Natural Gas Reservoirs, University of Michigan, 1955

[22] Pascal H, Quillian R G, Kingston J. Analysis of vertical fracture length and non-Darcy flow coefficient using variable rate tests. SPE 9348, 1980

[23] Koтяxoв Ф И. The Physical Foundation of Reservoir. Beijing: Petroleum Industry Press, 1958

[24] Li C. The research of natural gas small leakage & dispersion rule. Dissertation for the Doctoral Degree. Beijing: China University of Petroleum, 2011

[25] Saboorian-Jooybari H, Pourafshary P. Significance of non-Darcy flow effect in fractured tight reservoirs. Journal of Natural Gas Science and Engineering, 2015, 24: 132–143

[26] Wen Z, Huang G, Zhan H. An analytical solution for non-Darcian flow in a confined aquifer using the power law function. Advances in Water Resources, 2008, 31(1): 44–55

[27] Zhang Z Y, Nemcik J. Fluid flow regimes and nonlinear flow characteristics in deformable rock fractures. Journal of Hydrology (Amsterdam), 2013, 477: 139–151

[28] Brinkman H C. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Applied Scientific Research, 1949, 1: 27–34

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature