A novel POD reduced-order model based on EDFM for steady-state and transient heat transfer in fractured geothermal reservoir

International Journal of Heat and Mass Transfer 146 (2020) 118783

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

A novel POD reduced-order model based on EDFM for steady-state and transient heat transfer in fractured geothermal reservoir Tingyu Li a, Yanqing Gao b, Dongxu Han b,⇑, Fusheng Yang a,⇑, Bo Yu b a

School of Chemical Engineering and Technology, Xi’an Jiaotong University, Xi’an 710049, China School of Mechanical Engineering, Beijing Key Laboratory of Pipeline Critical Technology and Equipment for Deepwater Oil & Gas Development, Beijing Institute of Petrochemical Technology, Beijing 102617, China b

a r t i c l e

i n f o

Article history: Received 15 May 2019 Received in revised form 25 August 2019 Accepted 24 September 2019

Keywords: POD Reduced-order model Embedded discrete fracture model Porous media

a b s t r a c t A new proper orthogonal decomposition reduced-order model (POD-ROM) for fluid flow and heat transfer processes in a fractured geothermal reservoir is proposed. First, the primary governing equations are established based on the embedded discrete fracture model. Second, the POD-ROM is developed by projecting the discrete matrix equations instead of the primary governing equations onto the lowdimensional space spanned by the pressure and temperature basis functions. Compared with the typical POD modeling method, the new model reduction technique can easily deal with mass and energy transfer between the matrix and fractures, thus ensuring the mass and energy conservation of the whole system in the projection process. The new model does not require complex mathematical deduction and is easy to implement in practical applications. The accuracy and robustness of the model were carefully studied for three complex heat transfer cases. Calculation results indicated that the computational efficiency can be improved by at least 10–15 times while high reconstruction and prediction accuracy is maintained. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction With the increasing emphasis on environmental protection, the development and utilization of renewable energy are becoming increasingly important. Geothermal energy, as one type of renewable energy, has wide exploitation prospects and value owing to its high quality of thermal storage. Geothermal energy is generally extracted from great depths, and the permeability of the thermal storage region is relatively low because it contains few small natural fractures. It is often necessary to form an artificial fracture network via hydraulic fracturing to enhance the flow capacity of the working fluid. The fracture system provides the major flow path and heat exchange location. Therefore, the fractures must be explicitly characterized in the process of model development. Many fracture models, such as the discrete fracture network model (DFN) [1,2], discrete fracture model (DFM) [3,4], embedded discrete fracture model (EDFM) [5–9], have been proposed in the past few decades. Among them, the EDFM, which employs non-conforming grids to discretize the fracture network, exhibits good application ⇑ Corresponding authors. E-mail addresses: [email protected] (D. Han), [email protected] (F. Yang). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118783 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

prospects. The efficiency of numerical calculation can be improved. However, when the realistic length scales and complex fracture distributions are considered in the model domain, the computational cost is high. At present, researchers mainly use forward simulation to study geothermal extraction processes by establishing reliable mathematical models [1,4,10–15]. Although current computing power can solve this complicated problem, high-fidelity simulation is often computationally expensive. Therefore, it is necessary to establish an efficient method to accelerate the simulation of thermal extraction. Model reduction methods are effective for reducing the computational demand and have been successfully applied in many fields, e.g., natural convection [16], viscoelastic turbulent flows [17], crude oil pipelines [18], two-phase flows [19], etc. As indicated by the literature, the objective of model reduction methods is to construct a reduced-order model (ROM) in a low-dimensional space spanned by the basis functions for fast calculation. Proper orthogonal decomposition (POD) is a frequently used method to obtain the basis functions from high-fidelity simulation data. In recent years, POD-ROMs have been introduced into the field of underground porous media. We briefly review the representative works as follows.

2

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

Vermeulen et al. [20] applied the POD to establish two different ROMs for model inversion of an unsteady nonlinear groundwater flow. These two models were demonstrated to be computationally robust and more efficient than the classical inverse methodologies for several test cases. Siade et al. [21] investigated the effect of snapshot selection on the precision of a POD-ROM for groundwater flow. They proposed a simple exponential function to select snapshots at the appropriate times, which made the ROM significantly more accurate. Winton et al. [22] developed a new POD-ROM for saturated groundwater flow, which used sensitivities as the basis functions. Compared with the original model, the computational time was reduced by at least an order of magnitude. Boyce et al. [23] applied model reduction techniques to an unconfined groundwater-flow model. A new model reduction methodology based on Galerkin projection and the Newton formulation was proposed. Compared to the standard projection-based method, Ushijima et al. [24] proposed a fast algorithm for constructing the system matrices of a POD-ROM for a groundwater model with a high computational efficiency. Chaturantabut and Sorensen [25] developed a nonlinear reduced-order system for miscible viscous fingering in porous media by using the discrete empirical interpolation method (DEIM) in conjunction with POD. Numerical results indicated that the computational time for the POD-DEIM reduced systems was significantly reduced. Ghommem et al. [26] established an ROM for highly heterogeneous porous media by using dynamic mode decomposition (DMD) and POD. In different numerical examples, the DMD- and POD-based model exhibited a good capability to reproduce and predict the flow field. Rizzo et al. [27] proposed an adaptive POD-ROM for the advection-dispersion equation for porous media. The snapshot splitting technique (SST) was introduced to damp the temporal increase of the modeling error. As indicated by the foregoing review, many efforts have been directed toward reducing the amount of computational resources for the simulation of groundwater flows and other related flows in non-fractured porous media. However, as mentioned previously, fractures play a significant role in the subsurface flow and heat transfer processes in fractured geothermal reservoirs. Additionally, the scale and internal physical properties of fractures differ significantly from those of the surrounding matrix, increasing the computational challenges [28]. Therefore, it is necessary to consider the dimensionality reduction of the fracture system. Researchers have attempted to reduce the computational barrier for fractured porous media from different perspectives, mainly via multiscale models [29–32], and other technologies, such as machine learning [33] and regression-based ROMs [34,35]. According to the multiscale theory, Praditia et al. [29] established a multiscale formulation for coupled flow-heat equations for fractured geothermal reservoirs. Pressure and temperature multiscale basis functions were constructed by solving the local elliptic problems. Then, the nonlinear transient multiscale coarse-scale system was obtained. Compared with the fine-scale solutions, the developed multiscale method exhibited high accuracy. Vasilyeva et al. [30] proposed a multiscale method for simulating geothermal fractured reservoirs by using the generalized multiscale finite element method (GMsFEM). The multiscale basis functions are constructed by solving local problems in overlapping coarse regions. Accurate pressure and temperature solutions can be obtained by using a few multiscale basis functions. Nissen et al. [31] proposed a heterogeneity-preserving upscaling methodology for fractured geothermal reservoirs. The pressure equation is solved on finescale grids, which can obtain the conservative flux. The heat transport equation is discretized and solved on coarse grids. Numerical results indicated that the upscaling method produced accurate results compared with fine-scale solutions. On the basis of the nonlocal multicontinuum method (NLMC), Vasilyeva et al. [32]

presented an upscaling methodology for the single-phase flow and heat transport in fractured geothermal reservoirs. By constructing the local pressure and temperature basis functions in the oversampled domain, the required upscaled coarse-grid approximation can be obtained. The new proposed model can be viewed as an accurate general multiscale framework. Although multiscale methods have been applied to the dimensionality reduction of geothermal fractured systems, the topic has not been fully explored. Additionally, compared with the global model reduction techniques, the local model order reduction method is always more complex to implement, and it is sometimes difficult to control the numerical accuracy in practical engineering applications [36]. To the best knowledge of the authors, it is hard to see any work focusing on the POD-ROM for fracture-dominated flow and heat transfer problems. This is mainly because the commonly used POD Galerkin projection method cannot be directly applied to the coupled governing equations of matrix and fracture media, as it violates the conservation of mass and energy, and causes the failure of the POD modeling [37]. To address this challenge, in the present study, on the basis of embedded discrete fracture model (EDFM), we developed a PODROM for the fluid flow and heat transfer processes in a fractured geothermal reservoir for the first time. The model is suitable for steady-state and transient heat transfer conditions. The discrete matrix equations rather than primary governing equations are projected onto the low-dimensional space formed by a truncated series of known basis functions, which can avoid the non-zero mass and energy transfer between the matrix and the fracture media. The remainder of this paper is organized as follows. The primary governing equations, including the mass and energy conservation equations, are introduced in Section 2. The development of POD-ROMs for pressure and temperature fields is discussed in Section 3. The verification and application of the proposed models are presented in Section 4. Concluding remarks are presented in Section 5. 2. Primary governing equations 2.1. Flow equations The matrix and fractures are considered as two separate media in the EDFM, and they are governed by mass conservation equations expressed by Eqs. (1)–(10). The mass conservation equation for the matrix is

/m ct

X frim @pm Q þ Qm þ r  ðv m Þ ¼ W i @t

ð1Þ

The mass conservation equation for the ith fracture is

/fri ct

  @pfri þ r  v fri ¼ Q mfri þ Q frifrj þ Q fri W @t

ð2Þ

The velocity in the matrix and fractures is described by Darcy’s law.

vi ¼ 

i k 

lf



rpi  qf g ;

i ¼ m or fr

ð3Þ

The exchange flow Q frim between the matrix and the ith fracture is given as

Q frim ¼ CIkt ðpfri  pm Þ

ð4Þ

To ensure the conservation of the total flux between the matrix and the fractures, the exchange flow Q mfri can be defined as [38]

Z

Z Q mfri dA ¼  A

Q frim dV V

ð5Þ

3

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

The exchange flow Q frifrj between different fractures is given as

Q frifrj ¼ TIgt ðpfri  pfrj Þ

ð6Þ

Here, the superscripts m and fr represent the matrix and the fractures, respectively; ct is the coefficient of compressibility; / represents the porosity; p represents the pressure; t represents the time; QW is the source/sink term, which can be calculated using the Paceaman formulation [39–41]; k represents the permeability; lf represents the fluid viscosity; qf represents the fluid density; and g represents the gravitational acceleration. kt and gt represent the mean total mobility of the fluid at the matrix-fracture and fracture-fracture interfaces, respectively, and can be calculated as follows: m fri

kt;ij;r ¼

2k k  ij r  lf kmij þ kfri r

gt;k;r ¼

2k k  k r  frj lf kfri k þ kr

ð7Þ

fri frj

ð8Þ

Aij;k hdiij;k

ð9Þ

where Aij;k represents of fracture element k; hdiij;k represents the average normal distance from the matrix to the fracture, which is R defined as hdiij;k ¼ xk ðx0 Þdx0 =V ij ; xk represents the distance from the fracture within matrix element ij and V ij represents the volume of matrix element ij. TI represents the transmissivity between different fracture intersection elements, which can be calculated as [41,42]

TIi;j ¼

Ii  Ij ; Ii þ Ij

fr fr

Ii ¼ 2bi ki =



lf L i



ð10Þ

Z

Z

Emfri dA ¼  A

Efrim dV

ð14Þ

V

The energy exchange Efrifrj between different fractures is given as

      fri frj Efrifrj ¼ H Vfrifrj Vfrifrj h þ H Vfrifrj Vfrifrj h þ TInðT fri  T frj Þ ð15Þ In Eqs. (11)–(15), T represents the temperature; cpf represents the specific heat capacity of fluid; QE represents the injected or produced energy, which is calculated as Q E ¼ qf Q W cpf T; h represents

the specific fluid enthalpy, which is defined as h ¼ qf cpf T; V represents the relative velocity between the matrix and the fractures or between different fractures; and H(∙) represents the Heaviside step  i function. qcp eff is the effective physical parameter, which can be defined as



where the subscript ij indicates the matrix cell, and the subscript k and r represent the fracture segments in the ith and jth fractures, respectively. For discrete fracture element k overlapping with matrix element ij, the connectivity index CI is defined as [8]

CIij;k ¼

To ensure the conservation of the total flux between the matrix and the fractures, the exchange flow can be defined as

qcp

i eff

 i  i ¼ /i qf cpf þ ð1  /Þ qr cpr ;

with i ¼ m or fri

ð16Þ

i

keff represents the effective thermal conductivity, which can be calculated as i

i

i

keff ¼ /i kf þ ð1  /Þkr ;

with i ¼ m or fri

ð17Þ

where the subscripts f and r represent fluid and rock, respectively. f and n represent the mean heat conductivity at the matrixfracture and fracture-fracture interfaces, respectively, calculated via the harmonic average method.

fij;k ¼

nk;r ¼

fri 2km ij  kk

ð18Þ

fri km ij þ kk frj 2kfri k  kr

ð19Þ

frj kfri k þ kr

Here the subscript ij represents the matrix cell, and k and r denote the fracture segments.

fr

where bi represents the fracture aperture; Li represents the length of the fracture segment. 2.2. Energy equations Similarly, two sets of energy conservation equations (Eqs. (11) and (12)) are established to describe the matrix and fracture media, respectively. The energy exchange between them is described by Eqs. (13)–(15). The energy conservation equation of the matrix is



qcp

  X m @T m m m m m q c v  r T ¼ r  k r T Efrim þ Q m þ þ pf f eff E eff @t i

ð11Þ

The energy conservation equation of the ith fracture is



qcp

  fri @T fri fri þ Emfri þ Efrifrj þ Q fri þ qf cpf v fri  rT fri ¼ r  kfri E eff rT eff @t ð12Þ

The energy exchange Efrim between the matrix and the ith fractures is given as

    fri m Efrim ¼ H Vfrim Vfrim h þ H Vfrim Vfrim h þ CIfðT fri  T m Þ

ð13Þ

2.3. Discretization and solution of equations The matrix and fracture grids are presented in Fig. 1. In the EDFM, the fractures are reduced to one-dimensional (1D) line segments, which are embedded into the two-dimensional (2D) structured matrix grids. To achieve a high-fidelity numerical simulation, the finite volume method (FVM) is adopted to discretize the flow and energy equations. In the flow equations, discretization of the pressure gradient terms is implemented via the classical two-point flux approximation (TPFA) scheme. In the energy equations, the convection terms are discretized via the quadratic upwind interpolation of convective kinematics (QUICK) scheme. The second-order central differencing (CD) scheme is adopted for the diffusion terms. The fully implicit scheme is used for time discretization. For the matrix pressure equation, we obtain the following discrete form:

ct /m DDxDt y þ

   t   tþ1 tþ1 Dykw m þ Dxt pm pm  pm i;j i;j i;j  pi1;j

 tþ1  tþ1 Dxkn m m pm þ Dyt pm i;j  pi;j1 i;j  pi;jþ1   tþ1  tþ1 P ¼ CIk kt;i;j;k pfrk  pm þ J p  pm W i;j i;j

Dyket Dx



m pm i;j  piþ1;j

tþ1

þ

Dxkst Dy

Xi;j \Xk

ð20Þ

4

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

For the energy equation of fractures, the discrete form is given   k1=2  kþ1=2 tþ1  t   tþ1   þ qf cpf v T qcp eff DAft T frk  T frk bf  qf cpf v T bf tþ1 2b kkþ1=2  tþ1 k1=2  2bf keff f T fr  T frk1 þ Lk þLeffkþ1 T frk  T frkþ1 ðLk þLk1 Þ k ð Þ  tþ1 P  CIi;j fi;j;k T frk  T m i;j þ

Xk \Xi;j

¼

P Xk \Xr

 tþ1 frk TIk;r nk;r T frr s  Tk þ

  i P h  frkfrr  frkfrr frk  frr V H V hk þ H V frkfrr V frkfrr hr tþ1

Xk \Xr

ð23Þ

where the superscripts e, w, n and s represent the interfaces of the control volume. According to the discrete schemes used in Eqs. (20)–(23), the following two linear systems can be constructed. Fig. 1. Schematic of matrix and fracture grids for EDFM. Green point represents the production well, red points denote the matrix cells, and blue points are the fracture cells. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Amm Afm

Amf Aff

þ ¼

   t  tþ1 þ pfrk  pfrk kþ1=2

2bf kt

ðLk þLkþ1 Þ

k1=2

2bf kt

ðLk þLk1 Þ

 tþ1 pfrk  pfrk1

 tþ1 pfrk  pfrkþ1

 tþ1  tþ1 P TIk;r gt;k;r pfrs  pfrk  CIi;j kt;i;j;k pfrk  pm i;j

P Xk \Xr

Xk \Xi;j

ð21Þ Here, Dx and Dy represent the space intervals, and Dt represents the time interval; kt represents the mean total mobility at the control volume interface; bf represents the fracture aperture, and Lk represents the length of the kth fracture segment. Similarly, the discrete form for the energy equation of matrix can be expressed as



qcp

   t  tþ1 Tm  Tm i;j i;j



DxDy eff Dt

h e  w w þ qf cpf v ex T e Dy  qf cpf v w x T Dy itþ1  n  s þ qf cpf v ny T n Dx  qf cpf v sy T s Dx þ þ

Dykeeff Dx Dxkneff Dy

¼þ þ

fr

p

 ¼

; fr B

Bm

C mm

C mf

fm

ff

C

C

!

Tm T

fr



 ¼

Dm D



fr

ff

mm

ff

3. Construction of POD-ROMs for fluid flow and heat transfer As shown in Fig. 2, the modeling process of the POD-ROM comprises the following four steps. (1) Considering the 2D fluid flow and heat transfer problem in fractured porous media shown in Fig. 1, flow and energy equations are solved to obtain the representative sample data via the FVM. The sample assembles of the pressure and temperature fields are given in matrices Sp and ST. We do not need to consider samples at all time instances. Only representative pressure and temperature samples at the appropriate time-step, which capture the dominant information of the flow and heat transfer, should be put into the sampling matrices.

 tþ1 Dykw  tþ1 m m Tm þ Dxeff T m i;j  T iþ1;j i;j  T i1;j  tþ1 Dxks  tþ1 m m Tm þ Dyeff T m i;j  T i;jþ1 i;j  T i;j1

P Xi;j \Xk

  itþ1 P h  frm  frm fr m H V hk þ H V frm V frm hi;j V

Qm W c pf

ð22Þ

 tþ1 CIk fi;j;k T frk  T m i;j

Xi;j \Xk

þq



where A , A , C , and C represent the sub-matrix blocks, which contain the matrix-matrix and fracture-fracture transmissivities; Amf, Amf, Cmf and Cfm represent the off-diagonal sub-matrix blocks, which contain the matrix-fracture transmissivities; and Bm, Bf, Dm and Df are the terms on the right-hand side (RHS) of the discrete equations. The detailed definitions of the coefficient matrices are given in the Appendix A. The linear systems expressed by Eq. (24) can be solved via a direct method or an iterative method. For large matrix systems, this incurs a significant computational burden [43]. To resolve this difficulty, POD-ROMs are developed, as described in detail in Section 3.

For the fracture pressure equation, we have: A

pm

ð24Þ mm

ct /fr Dft

!

 tþ1 Tm i;j

ð25Þ Here, tn represents the nth sampling time-step; i and j represent the total numbers of control points in the matrix and fractures, respectively.

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

For convenience, the sampling matrices Sp and ST can be rewritten as follows.

Sa ¼ ½a1 ; a2 ;    ; an1 ; an rn ;

a ¼ p or T

ð26Þ

 Tr where the pressure sampling matrix is p ¼ pm ; pfr and the tem Tr m fr perature sampling matrix is T ¼ T ; T . The number of rows in the sampling matrix is r ¼ i þ j. The transpose operator is defined as Tr (to distinguish it from the temperature superscript T). (2) According to the sampling matrix Sa, the basis functions (or POD modes) can be obtained via the ‘‘snapshot” method [44–46] (if r > n) or the singular value decomposition method (SVD) (if r < n). Generally, the total number of control points is far greater than the number of samples; thus r > n. Therefore, the ‘‘snapshot” method is performed to obtain the eigenvalues and eigenvectors. As indicated by Eqs. (27) and (28), the inner-product matrix Ma is first constructed, and the eigenvalues and corresponding eigenvectors are then calculated.

Ma ¼ Sa Tr Sa

ð27Þ

M a V a ¼ ka V a

ð28Þ

Here ka represents the eigenvalues, and Va represents the corresponding eigenvectors. Finally, the pressure and temperature basis functions can be computed as follows.

/a ¼

Sa V a k Sa V a k

ð29Þ

 Tr fr where the basis functions are /a ¼ /m . a ; /a The use of more basis functions in the reconstruction processes does not mean that the accuracy of the calculation is higher. The number of selected basis functions is often determined via trial calculations. In this study, we use the following two criteria to help us truncate basis functions.

5

(a) Energy contribution of the ith POD basis function

b i ¼ ki =

n X

kk

ð30Þ

k¼1

(b) Accumulative energy contribution of the first m POD basis functions

u¼

m X

kk =

k¼1

n X

kk ;

m6n

ð31Þ

k¼1

The magnitude of the energy contribution helps us identify which basis functions capture the key information of the original physical problem. When the accumulative energy contribution is approximately 100%, the number of POD basis functions can be truncated. (3) According to the POD theory, the pressure and temperature fields can be expressed as a linear combination of basis functions and spectral coefficients.

p¼

m X

apk /pk ;

ð32aÞ

aTk /Tk

ð32bÞ

k¼1

T¼

w X k¼1

where ap and aT represent the spectral coefficients of the pressure and temperature fields, respectively, and m and w represent the number of basis functions selected. For the sake of simplicity, Eq. (24) can be rewritten in a more concise form.

Ann pn1 ¼ Bn1

ð33aÞ

Cnn Tn1 ¼ Dn1

ð33bÞ

Here A and C represent the discrete coefficient matrices, and B and D represent the RHS terms of the discrete matrix equations. As shown in Fig. 2, in the traditional modeling process of the POD-ROM, the primary governing equations are directly projected

Fig. 2. Modeling process of POD-ROMs.

6

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

Fig. 3. The traditional and proposed modeling method of ROM.

y

y

o

x

o (a)

x (b)

y

o

x (c)

Fig. 4. Computational domains and boundary conditions: a for steady-state case, b and c for transient cases.

7

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

onto the low-dimensional space formed by the basis functions [16–18]. After complex equation deductions, the ROM can be obtained. However, this modeling method is not suitable for the fracture-dominated flow and heat transfer problem expressed by Eqs. (1)–(2) and (11)–(12). Specifically, there are two ways (procedures (a) and (b) shown in Fig. 3) to build a ROM, if the traditional POD modeling process is implemented [37]. Procedure (a): The governing equations of the matrix and fracture media are projected onto their respective POD basis functions to construct the corresponding ROM. The final POD-ROM is obtained by adding them together. However, this approach introduces additional mass and energy transfer terms into the PODROM, leading to modeling failure. Procedure (b): Before modeling, the governing equations of matrix and fracture media are added to eliminate the mass and energy coupling terms artificially. Then, the composite equations can be projected onto the POD basis functions related only to the matrix or fracture media. Testing is needed to determine which ROM has higher accuracy. This approach is not universal, as it severely influenced by the artificial and empirical factor. To overcome these hurdles, a new model reduction strategy (procedure (c), as shown in Fig. 3) is proposed. By projecting the discrete matrix equations instead of the primary governing equations onto the low-dimensional space spanned by the first m pressure and w temperature basis functions, the model reduction equations can be established.

 

/pnm /Tnw

Tr

 Tr Ann /pnm apm1 ¼ /pnm Bn1

ð34aÞ

Tr

 Tr Cnn /Tnw aTw1 ¼ /Tnw Dn1

ð34bÞ

case (Fig. 4(a)) is designed for steady-state flow and heat transfer conditions. Fig. 4(b) and (c) are applied for transient cases to simulate the thermal extraction processes of fractured geothermal reservoirs. The model size is 500 m  500 m. The other main model parameters are presented in Table 1. 4.1. Steady-state heat transfer case with 8 discrete fractures As shown in Fig. 4(a), the square matrix contains eight discrete fractures. The grid number of the matrix is 10,000, and the Table 1 Main computational parameters of model. Model Parameters

Values

Units

Permeability of matrix and fracture Density of matrix and fluid Porosity of matrix and fracture Thermal conductivity of matrix and fluid Thermal capacity of matrix and fluid Fluid viscosity Fracture aperture

1015, 109 2600, 1000 0.2, 0.3 2.0, 0.5 1000, 4000 103 103

m2 kg/m3 – W/(mK) J/(kgK) Pas m

Therefore, the POD-ROMs of the pressure and temperature fields can be directly constructed. 







Amm apm1 ¼ Bm1

ð35aÞ

Cww aTw1 ¼ Dw1 

ð35bÞ 

where Amm and Cww represent the coefficient matrices of the    Tr  Tr ROMs (Amm ¼ /pnm Ann /pnm , Cww ¼ /Tnw Cnn /Tnw ), and 



Bm1 and Dw1 represent the RHS terms of the matrix equations    Tr  Tr (Bm1 ¼ /pnm Bn1 ,Dw1 ¼ /Tnw Dn1 ). Clearly, Eq. (35) is a linear matrix equation, which can be directly solved to obtain the spectral coefficients by using the backslash operator in MATLAB. Compared with Eq. (33), the number of degrees of freedom of Eq. (34) is reduced. The simulation efficiency can be significantly increased.

(a) Pressure

(4) The reconstruction processes of approximate solutions mainly include two steps. The first step is to use Eq. (34a) to compute the spectral coefficient of pressure (apk ). Using Eq. (32a), the pressure field is first reconstructed. Accordingly, Darcy’s velocity can be computed, which is used to calculate the matrix Cnn . The second step is to solve Eq. (35b) to obtain the spectral coefficient of temperature (aTk ). Subsequently, the approximate temperature field can be reconstructed using Eq. (32b). 4. Verification and application Even though Eq. (35) is deduced from transient governing equations, it is still suitable for steady-state conditions. Therefore, to verify the accuracy and robustness of the proposed POD-ROM, we design three synthetic cases, as shown in Fig. 4. The first test

(b) Temperature Fig. 5. Energy contribution and cumulative energy contribution for pressure and temperature basis functions.

8

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

fractures are discretized by 652 grids. The pressure and temperature boundary conditions on the four sides of the matrix are applied with the following nonlinear functions in the x- and ydirections [47].

Left side : pl ¼ a1 þ a2 siny;

T l ¼ a1 þ a2 cosy

Right side : pr ¼ a3 þ a4 siny; Topside : pt ¼ a5 þ a6 sinx;

ð36Þ

T r ¼ a3 þ a4 cosy

ð37Þ

T t ¼ a5 þ a6 cosx

Bottomside : pb ¼ a7 þ a8 sinx;

ð38Þ

T b ¼ a7 þ a8 cosx

ð39Þ

There are eight variables, i.e., a1–a8, in Eqs. (36)–(39). To obtain a series of pressure and temperature samples, a1–a8 are selected with two values; 0 and 106 for pressure and 0 and 100 for temperature. Therefore, the total numbers of samples is 28 = 256 for both the pressure and temperature. The basis functions of the pressure and temperature are then obtained via the POD method. According to Eqs. (30) and (31), the energy contribution and cumulative energy contribution for the pressure and temperature basis functions are presented in Fig. 5. For the pressure basis functions, the first seven basis functions are more energetic than the others. The total energy contribution of the first eight basis functions is very close to 100%. Therefore, the first eight basis functions are chosen for pressure reconstruction. For the temperature

(a) φ1p (89.5%)

(d) φ1T (92.3%)

basis functions, the energy contribution of the first 10 basis function is significantly higher than of the others. The cumulative energy contribution of the first 14 basis functions is 99.99%. Thus, the first 14 basis functions are selected for temperature reconstruction. To show the basis function more intuitively, the spatial distributions of different pressure and temperature basis functions are displayed in Fig. 6. High-energy basis functions often capture the important information affecting flow and heat transfer processes, such as boundary conditions and fractures, as shown in /p1 , /p3 , /T1 and /T3 . The low-energy basis functions depicted in /p8 and /T8 mainly capture the information of the boundary area where the pressure and temperature change sharply. To quantitatively verify the reconstruction and prediction accuracy of the developed POD-ROM, we define the average relative deviation between the POD and the FVM.

eF

 PNtot 

POD

 F FVM

l l¼1 F l



; ¼ PNtot FVM

l¼1 F l

F ¼ p or T

ð40Þ

Here, Ntot represents the total number of control points, and the superscripts POD and FVM represent the results of POD-ROM and the full-order model, respectively.

(b) φ3p (3.8%)

(c) φ8p (0.053%)

(e) φ3T (2.5%)

(f) φ8T (0.064%)

Fig. 6. Space distributions and energy contributions of different pressure and temperature basis functions: the first row represents the pressure basis functions, and the second row represents the temperature basis function.

Table 2 Deviations between the POD and the FVM for reconstructions of original pressure and temperature samples. Presssure

pl (MPa)

pr (MPa)

pt (MPa)

pb (MPa)

Deviations

Sample Sample Sample Sample

siny siny 0 1

siny 1 + siny siny 0

sinx 1 + sinx 0 0

sinx 1 + sinx sinx 1

1.40  109% 2.10  109% 1.51  109% 1.39  109%

Tl (100 °C) cosy cosy 0 1

Tr (100 °C) cosy 1 + cosy cosy 0

Tt (100 °C) cosx 1 + cosx 0 0

Tb (100 °C) cosx 1 + cosx cosx 1

1.13  105% 3.56  105% 1.67  105% 2.33  105%

1 2 3 4

Temperature Sample 1 Sample 2 Sample 3 Sample 4

9

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

First, the reconstruction accuracy of the POD-ROM for the original pressure and temperature samples is verified. As indicated by Table 2, eight original samples with corresponding pressure and temperature boundary conditions are randomly selected. Among them, we clearly see that the maximum deviations for the pressure and temperature reconstructions are 2.10  109% and 3.56  105%, respectively, which are acceptable for practical calculations. Fig. 7 shows a comparison of the reconstruction results for the original pressure and temperature Sample 2. The reconstructed pressure and temperature of the POD-ROM agree well with the results of FVM in both the boundary and interior regions. Therefore, the reconstruction accuracy of the ROM is very high.

Second, the prediction accuracy of the ROM for out of original pressure and temperature samples is verified. We design eight prediction conditions, as shown in Table 3. The boundary conditions of the pressure and temperature used in Prediction 1 are the same, whereas the other three prediction cases have different boundary conditions. The left and top conditions used in Prediction 3 are within the scope of the original samples, and the remaining conditions are beyond the scope. Among these prediction cases, the maximum deviation between the POD and the FVM is 1.12  108% for the pressure, and 5.46  104% for the temperature, indicating the applicability of the proposed POD-ROM. Fig. 8 shows two sets of results for pressure and temperature fields.

(a) Pressure

(b) Temperature Fig. 7. Comparison of reconstruction results for original pressure and temperature Sample 2. Left: FVM, Right: POD.

Table 3 Deviations between the POD and the FVM for predictions of pressure and temperature. Pressure Prediction Prediction Prediction Prediction

1 2 3 4

Temperature Prediction 1 Prediction 2 Prediction 3 Prediction 4

pl (MPa)

pr (MPa)

pt (MPa)

pb (MPa)

Deviations

2.7 3 + 2siny 0 2.5siny

1.6 2 + 3siny 3 2.7

2.9 2.5 + 3.5sinx 0 2.4 + 2.3sinx

1.4 3.7 + 2.6sinx 3sinx 1.5 + 1.6sinx

5.59  109% 1.12  108% 5.77  109% 4.66  109%

Tl (100 °C) 2.7 2 + 3cosy 0 2.7cosy

Tr (100 °C) 1.6 1.5cosy 1.7 2.5

Tt (100 °C) 2.9 2.5 + 3.5cosx 0 2.3

Tb (100 °C) 1.4 3.7 + 2.6cosx 0.5sinx 1.6cosx

4.48  104% 5.46  104% 1.05  104% 1.78  104%

10

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

(a) Pressure

(b) Temperature Fig. 8. Comparisons of prediction results of pressure and temperature fields. Left: FVM, Right: POD.

(a) φ1T (98.9%)

(d) φ4T (0.011%)

(b) φ2T (0.95%)

(e) φ5T (1.5×10-3%)

(c) φ3T (0.089%)

(f ) φ6T (1.4×10-4%)

Fig. 9. Spatial distribution and energy contribution of the first six temperature basis functions.

11

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

The results of the FVM and the POD are shown on the left, and right sides, respectively. These two groups of results exhibit almost identical pressure and temperature distributions. Therefore, the proposed POD-ROM is effective for reproducing and predicting the flow and heat transfer processes in this steady-state case. 4.2. Transient thermal extraction case 1 with 7 discrete fractures In this subsection, we mainly discuss the applicability of the proposed POD-ROM for the transient thermal extraction process. As shown in Fig. 4(b), the computational domain contains seven discrete fractures, which characterizes the fracture network after hydraulic fracturing. The grid numbers of the matrix and fracture are 10,000 and 306, respectively. The top and bottom sides are both subjected to an impermeable boundary condition for pressure. The adiabatic boundary conditions for temperature are applied to the top, bottom and right boundaries. The left and right boundaries simulate the injection well and production well, respectively. The injection pressure and temperature are 30 MPa and 20 °C, respectively. The right boundary simulates the production well with a production pressure of 1 MPa. The initial reservoir temperature is 180 °C. The time-step is 0.1 years, and the total simulation time is 30 years. In this case, the pressure field remains almost stable after approximately 10 d of injection. The temperature field changes slowly; therefore, we focus on the reconstruction and prediction accuracy of the POD-ROM for the transient temperature field. One straightforward way to obtain a representative snapshot containing dominant characteristics of the physical problems is to apply various boundary conditions and parameter settings. For this case, to capture the unsteady process of temperature, the sampling begins at t = t0, and the sampling interval is set as one year; thus, a single temperature sample is selected every other year. The sampling method for the pressure field is similar. Therefore, the total number of pressure and temperature samples is 60. The spatial distribution and energy contribution of the temperature basis functions after sample decomposition using the "snapshot"

(a) 5 years

(b) 20 years

method are presented in Fig. 9. The first basis function occupies 98.9% of the total energy. The cumulative energy contribution of the first six basis functions is approximately 100%. Therefore, the first six basis functions are provided to the POD-ROM for temperature reconstructions. The first two pressure basis functions are selected to calculate the pressure. Fig. 10 presents the temperature field reconstructed using the POD method at 5, 20, and 30 years, in comparison with the FVM results. We observe no differences between the two groups of results. According to Eq. (40), the maximum relative deviation among the three reconstruction results is 3.43  102%, which can be neglected in practical applications. To further test the predictive ability of the proposed POD-ROM, we compare the temperature distribution of the right boundary with that for the FVM after 40 years of thermal extraction, as illustrated in Fig. 11. The maxi-

Fig. 11. Comparison of predicted temperature field between the FVM and the POD at 40 years.

(c) 30 years

Fig. 10. Comparisons of reconstructed temperature field at various years: the first row represents the results of FVM, and the second row represents the results of POD.

12

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

mum relative deviation between the FVM and the POD-ROM is 0.11%, indicating the good predictive ability of ROM. 4.3. Transient thermal extraction case 2 with 520 discrete fractures Here, we investigate the applicability of the POD-ROM for the transient thermal extraction process with wells. As shown in

Fig. 4(c), the computational domain contains injection and production wells and 520 discrete fractures. The coordinates of the injection well are (105 m, 105 m), and the production well is located at (395 m, 395 m). The distance between the wells is approximately 410 m, and the wellbore radius is 0.1 m. The numerical calculation parameters used in this case refer to the Yangbajing geothermal field [48]. The injection flow rate is 20 kg/

Fig. 12. Spatial distribution and energy contribution of the first three pressure and temperature basis functions.

(a) 5 years

(b) 20 years

(c) 30 years

Fig. 13. Comparisons of reconstructed temperature field at various years: the first row represents the results of FVM, and the second row represents the results of POD.

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

s, and the production pressure is 5 MPa. The injection temperature is 60 °C, and the initial reservoir temperature is 240 °C. The impermeable boundary conditions are applied on the four sides. The adiabatic boundary conditions are applied around the overall model. The grid numbers of matrix and fracture are 10,000 and 14,099, respectively. The total number of computational grids is 24,099. The time-step is 0.1 years, and we perform transient thermal extraction for a period of 30 years. The sampling process is similar to that used in transient case 2. The total number of pressure and temperature samples is 60. The spatial distribution and energy contribution of the first three pressure and temperature basis functions are presented in Fig. 12. The first eight basis functions capture almost 100% of the transient heat transfer process. Therefore, the ninth to the last basis function can be neglected in the temperature reconstruction process. The first five pressure basis functions are provided to reconstruct the seepage field. For the prediction process, the number of basis functions selected needs to be examined in practice. Compared with the original temperature samples, the reconstructed transient thermal extraction processes obtained using the POD-ROM are presented in Fig. 13. Clearly, the proposed model perfectly reproduces the whole transient heat transfer process using the eight aforementioned temperature basis functions. The maximum relative deviation between the FVM and the POD-ROM is 8.70  102%. A good POD-ROM can not only reproduce sample data, but also predict unknown information in a certain range. Fig. 14 presents the reconstructed and predicted outlet temperature after 40 years of thermal extraction. When the eight basis functions are used, the ROM cannot only reconstruct the original sample data with high accuracy but also predict the outlet temperature for the next 10 years. The prediction error is 0.83%. The computational efficiency of the POD-ROM is of great concern to us. For the steady-state case, the full-order model must solve 21,304 equations simultaneously, whereas the ROM only needs to solve 22 equations. Thus, the speed of numerical calculation can be increased by at least 10 times. For the transient case1, the number of degrees of freedom for the full-order model and ROM are 20,612 and 8, respectively. For the transient case2, compared with the 48,198 equations of the original model, the ROM only requires the solution of 13 equations. The computational efficiency can be increased by at least 15 times. Therefore, the proposed model reduction approach significantly reduces the computational effort.

13

5. Conclusions We proposed a POD-ROM based on the EDFM for steady-state and transient flow and heat transfer processes in fractured geothermal reservoirs. The main conclusions are as follows. (1) By projecting the discrete matrix equations instead of the primary governing equations onto the low-dimensional space formed by the pressure and temperature basis functions, the established POD-ROM can easily handle the coupling terms between the matrix and fracture media, and successfully avoid the non-zero mass and energy transfer between them, in contrast to the traditional POD projection method. (2) The developed POD-ROM is applied to reconstruct and predict the flow and heat transfer processes in three artificial fractured geothermal reservoirs. The POD-ROM cannot only reproduce the original sample data perfectly, but also predict the flow and heat transfer processes in fractured porous media within a certain range. For the steady-state case, the reconstruction and prediction accuracies of the POD-ROM are 2.10  109% and 1.12  108%, respectively, for pressure field and 3.56  105% and 5.46  104%, respectively, for the temperature field. Among the transient cases, the minimum reconstruction accuracy for the whole temperature field is 8.7  102%, and the prediction accuracy for the outlet temperature is 0.83%. The prediction accuracy of the model for steady-state heat transfer conditions is higher than that for transient conditions. In both cases, the accuracy of the model is acceptable for real applications. (3) The proposed POD-ROM has not only a high reconstruction and prediction accuracy but also a high calculation speed. Compared with the original full-order model solved by using the FVM, the computational efficiency can be increased by a factor of at least 10–15 using the ROM. The results indicate that POD combined with Galerkin projection is an efficient model reduction technique, which can significantly reduce the computational effort for fractured porous media. Although the proposed ROM is validated by 2D examples, the new model reduction method can be completely applied to the fast simulation of the three-dimensional (3D) thermal extraction process in a real geothermal reservoir, which should be investigated in a future study. Predictably, the acceleration effect will be more obvious under 3D conditions. The new proposed modeling method is general, and is not limited to a particular number of coupled equations. However, the POD method has limitations: the acquisition of basis functions needs to be based on a large amount of sample data, and the number of optimal basis functions must be determine via trial and error. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgments

Reconstruction

Prediction

Fig. 14. Comparison of predicted outlet temperature with eight POD basis functions.

The study is supported by the National Natural Science Foundation of China (No.51936001, No. 51706021), the Project of Construction of Innovative Teams and Teacher Career Development for Universities and Colleges under Beijing Municipality (No. IDHT20170507), the Beijing Youth Talent Support Program (CIT&TCD201804037), and the Program of Great Wall Scholar (No. CIT&TCD20180313). The authors thank the reviewers for their constructive and professional comments.

14

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783

Appendix A

F FUD ¼

In this appendix, we present the detailed definitions of the coefficient matrices for Eq. (24) in Section 2.3. The matrix flow (Eq. (20)) can be rewritten as follows. m m m m m m m m m frm fr am pk þ b P pi;j ¼ aW pi1;j þ aE piþ1;j þ aS pi;j1 þ aN pi;jþ1 þ ac

m

where

afrm c

Dykw t Dx

Dyke

t ; am E ¼ Dx ; P ¼ CIk kt;i;j;k

Dxkst Dy

am S ¼

;

am N ¼

Dxknt Dy

;





qf cpf v x  qf cpf v x e am S

 i h   i Dy þ qf cpf v y  qf cpf v y Dx w

n

s

FUD m m am þ am þ afrm þ ad þ awell þ am c P ¼ aW þ aE þ N þF p0 ,  t m m m QUICK . b ¼ ap0 T i;j þ s

Similarly, for the fracture energy (Eq. (23)), we have:

ðA:1Þ

am W ¼

h

m DxDy am P0 ¼ ct / Dt ;

afrp T frk ¼ afrk1 T frk1 þ afrkþ1 T frkþ1 þ afrkfrr T frs þ amfr Tm i;j þ b c c afrk1 ¼

k1=2

2bf keff

ðLk þLk1 Þ

  þ max qf cpf v ; 0

k1=2

fr

ðA:4Þ

bf ,

  afrkþ1 ¼ Lk þLk1 þ max qf cpf v ; 0 bf , ð Þ kþ1=2 P P mfr frkfrr ac ¼ CIi;j fi;j;k , ac ¼ TIk;r nk;r , kþ1=2

2bf keff

Xi;j \Xk

mfr m m m m am þ am P ¼ aW þ aE þ aS þ aN þ ac p0 þ J  t m m m b ¼ ap0 pi;j þ JpW

Xk \Xi;j

afrd ¼

For the fracture flow (Eq. (21)), we have another discrete form.

afrp pfrk ¼ afrk1 pfrk1 þ afrkþ1 pfrkþ1 þ afrkfrr pfrs þ amfr pm c c i;j þ b

fr

ðA:2Þ

Xk \Xr

  afrp0 ¼ qcp eff afrp ¼ fr

b ¼

where

Xk \Xr

 P   frkfrr  frkfrr V H V qc p ,

DxDy , Dt fr fr ak1=2 þ akþ1=2 þ amfr þ c  t fr fr QUICK QUICK ap0 T k þ s .s

afrkfrr þ afrd þ afrp0 ,c represents the correction terms of

the QUICK formulation. The specific theory is presented in [49]. afrk1 ¼ afrkfrr c

2b k

k1=2

2b k

kþ1=2

f t ; afrkþ1 ¼ LkfþLt k1 ; ðLk þLk1 Þ ð Þ P P ¼ TIk;r gt;k;r ; amfr ¼ CIi;j kt;i;j;k c

Xk \Xr

According to Eqs. (A.3) and (A.4), the specific coefficient matrices can be obtained.

Xk \Xi;j

C mm C fm

afrp ¼ afrk1 þ afrkþ1 þ afrkfrr þ amfr þ ct /fr DDxDt y c c  t fr b ¼ ct /fr DDxDt y pfrk

C mf C ff

!

Tm T

fr



 ¼

Dm D



fr

where Eqs. (A.1) and (A.2) can be further constructed as the following discrete system.

Amm

Amf

fm

ff

A

!

pm



fr

p

A

 ¼

Bm



fr

B

where  m n o n o i ¼ j mf aP ; frm m Amm ¼ am ¼ ; ; A ¼ amf , a , amf m m m ij ij ij ij ¼ ac am ; a ;a ;a ; i–j W E S N ( n o n o afrp ; i¼j mfr Afm ¼ afm ; Aff ¼ affij , affij ¼ , afm ; ij ij ¼ ac afrk1 ; afrkþ1 ; i–j n o m m fr fr fr Bm ¼ bi1 , bi1 ¼ bm ; Bfr ¼ bi1 , bi1 ¼ b . For the matrix energy (Eq. (22)), we can obtain the following discrete system. m m m m m m m m m frm fr am Tk þ b P T i;j ¼ aW T i1;j þ aE T iþ1;j þ aS T i;j1 þ aN T i;jþ1 þ ac

m

ðA:3Þ where   þ max qf cpf v x ; 0 Dy, w   Dykeeff m aE ¼ Dx þ max qf cpf v x ; 0 Dy,   e Dxkseff am S ¼ Dy þ max qf c pf v x ; 0 Dx, s   Dxkneff ¼ þ max  q c v ; 0 Dx, am f pf x N Dy n   P P ¼ CIk fi;j;k , ad ¼ H V frm V frm qf cpf , afrm c

am W ¼

Dykw eff Dx

Xi;j \Xk

awell ¼ q

m Qm W cpf ,ap0



¼ qc p

Xi;j \Xk



eff

DxDy , Dt

 m n o n o i¼j aP ; m ; C mf ¼ cmf C mm ¼ cm ¼ , c , m m m ij ij ij am ; a ; a ; a ; i–j W E n o S N n o , cfm ¼ afrm ; C fm ¼ cfm ¼ amfr ; C ff ¼ cffij , cmf c c ij ij ij ( n o m m afrp ; i¼j fr cffij ¼ ; Dm ¼ di1 , di1 ¼ bm ; Dfr ¼ di1 , fr fr ak1 ; akþ1 ; i–j fr

fr

di1 ¼ b .

References [1] T. Doe, R. Mclaren, W. Dershowitz, Discrete fracture network simulations of enhanced geothermal systems, Proceedings, Thirty-Ninth Workshop on Geothermal Reservoir Engineering, Stanford University, Stanford, California, 2014. [2] J.D. Hyman, S. Karra, N. Makedonska, et al., dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport, Comput. Geosci. 84 (2015) 10–19. [3] T.T. Garipov, M. Karimi-Fard, H.A. Tchelepi, Discrete fracture model for coupled flow and geomechanics, Comput. Geosci. 20 (1) (2016) 149–160. [4] J. Yao, X. Zhang, Z.X. Sun, et al., Numerical simulation of the heat extraction in 3D-EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures model, Geothermics 74 (2018) 19–34. [5] S.H. Lee, M.F. Lough, C.L. Jensen, Hierarchical modeling of flow in naturally fractured formations with multiple length scales, Water Resour. Res. 37 (3) (2001) 443–455. [6] L. Li, S.H. Lee, Efficient field-scale simulation of black oil in naturally fractured reservoir through discrete fracture networks and homogenized media, SPE Reservoir Eval. Eng. 11 (4) (2008) 750–758. [7] S.H. Lee, C.L. Jensen, M.F. Lough, Efficient finite-difference model for flow in a reservoir with multiple length-scale fractures, SPE J. 3 (5) (2000) 268–275. [8] H. Hajibeygi, D. Karvounis, P. Jenny, A hierarchical fracture model for the iterative multiscale finite volume method, J. Comput. Phys. 230 (24) (2011) 8729–8743. [9] J. Gunnar, S.A. Miller, On the role of thermal stresses during hydraulic stimulation of geothermal reservoirs, Geofluids 2017 (2017) 1–15.

T. Li et al. / International Journal of Heat and Mass Transfer 146 (2020) 118783 [10] A.R. Shaik, S.S. Rahman, N.H. Tran, et al., Numerical simulation of fluid-rock coupling heat transfer in naturally fractured geothermal system, Appl. Therm. Eng. 31 (10) (2011) 1600–1606. [11] F.Z. Zhang, P.X. Jiang, R.N. Xu, System thermodynamic performance comparison of CO2-EGS and water-EGS systems, Appl. Therm. Eng. 61 (2013) 236–244. [12] Z.Q. Qu, W. Zhang, T.K. Guo, Influence of different fracture morphology on heat mining performance of enhanced geothermal systems based on COMSOL, Int. J. Hydrogen Energy 42 (2017) 18263–18278. [13] Y.C. Zeng, J.M. Zhan, N.G. Wu, et al., Numerical investigation of electricity generation potential from fractured granite reservoir by water circulating through three horizontal wells at Yangbajing geothermal field, Appl. Therm. Eng. 104 (2016) 1–15. [14] Y. Chen, G.W. Ma, H.D. Wang, et al., Evaluation of geothermal development in fractured hot dry rock based on three dimensional unified pipe-network method, Appl. Therm. Eng. 136 (2018) 219–228. [15] H.S. Vik, S. Salimzadeh, H.M. Nick, Heat recovery from multiple-fracture enhanced geothermal systems: the effect of thermoelastic fracture interactions, Renewable Energy 121 (2018) 606–622. [16] D.X. Han, B. Yu, J.J. Chen, et al., POD reduced-order model for steady natural convection based on a body-fitted coordinate, Int. Commun. Heat Mass Transfer 68 (2015) 104–113. [17] J.J. Chen, D.X. Han, B. Yu, et al., A POD-Galerkin reduced-order model for isotropic viscoelastic turbulent flow, Int. Commun. Heat Mass Transfer 84 (2017) 121–133. [18] B. Yu, G.J. Yu, Z.Z. Cao, et al., Fast calculation of the soil temperature field around a buried oil pipeline using a body-fitted coordinates-based PODgalerkin reduced-order model, Numer. Heat Transfer, Part A: Appl. 63 (10) (2013) 776–794. [19] M. Ghasemi, E. Gildin, Model order reduction in porous media flow simulation using quadratic bilinear formulation, Comput. Geosci. 20 (3) (2016) 723–735. [20] P.T.M. Vermeulen, C.B.M.T. Stroet, A.W. Heemink, Model inversion of transient nonlinear groundwater flow models using model reduction, Water Resour. Res. 42 (9) (2006) 350–359. [21] A.J. Siade, M. Putti, W.W.G. Yeh, Snapshot selection for groundwater model reduction using proper orthogonal decomposition, Water Resour. Res. 46 (8) (2010) 2657–2662. [22] C. Winton, J. Pettway, C.T. Kelley, et al., Application of Proper Orthogonal Decomposition (POD) to inverse problems in saturated groundwater flow, Adv. Water Resour. 34 (12) (2011) 1519–1526. [23] S.E. Boyce, T. Nishikawa, W.W.G. Yeh, Reduced order modeling of the Newton formulation of MODFLOW to solve unconfined groundwater flow, Adv. Water Resour. 83 (2015) 250–262. [24] T.T. Ushijima, W.W.G. Yeh, A proposed Fast algorithm to construct the system matrices for a reduced-order groundwater model, Adv. Water Resour. 102 (2017) 68–83. [25] S. Chaturantabut, D.C. Sorensen, Application of POD and DEIM to dimension reduction of nonlinear miscible viscous fingering in porous media, 2009. Available: https://hdl.handle.net/1911/102127. [26] M. Ghommem, M. Presho, V.M. Calo, et al., Mode decomposition methods for flows in high-contrast porous media. Global-local approach, J. Comput. Phys. 253 (2013) 226–238. [27] C.B. Rizzo, F.P.J. De Barros, S. Perotto, et al., Adaptive POD model reduction for solute transport in heterogeneous porous media, Comput. Geosci. 22 (2017) 297–308. [28] S. Shah, O. Møyner, M. Tene, et al., The multiscale restriction smoothed basis method for fractured porous media (F-MsRSB), J. Comput. Phys. 318 (2016) 36–57.

15

[29] T. Praditia, R. Helmig, H. Hajibeygi, Multiscale formulation for coupled flowheat equations arising from single-phase flow in fractured geothermal reservoirs, Comput. Geosci. 22 (5) (2018) 1305–1322. [30] M. Vasilyeva, M. Babaei, E.T. Chung, et al., Multiscale modeling of heat and mass transfer in fractured media for enhanced geothermal systems applications, Appl. Math. Model. 67 (2019) 159–178. [31] A. Nissen, E. Keilegavlen, T.H. Sandve, et al., Heterogeneity preserving upscaling for heat transport in fractured geothermal reservoirs, Comput. Geosci. 22 (2) (2018) 451–467. [32] M. Vasilyeva, M. Babaei, E.T. Chung, et al., Upscaling of the single-phase flow and heat transport in fractured geothermal reservoirs using nonlocal multicontinuum method, Comput. Geosci. 4 (2019) 1–15. [33] S. Srinivasan, S. Karra, J. Hyman, et al., Model reduction for fractured porous media: a machine learning approach for identifying main flow pathways, Comput. Geosci. 23 (2019) 617–629. [34] M.K. Mudunuru, S. Kelkar, S. Karra, et al., Reduced-order models to predict thermal output for enhanced geothermal systems, Proceedings of 41st Stanford Geothermal Workshop, Stanford University, Stanford, CA, USA, 2016. [35] M.K. Mudunuru, S. Karra, D.R. Harp, et al., Regression-based reduced-order models to predict transient thermal output for enhanced geothermal systems, Geothermics 70 (2017) 92–205. [36] J.F. Li, T. Zhang, S.Y. Sun, et al., Numerical investigation of the POD reducedorder model for fast predictions of two-phase flows in porous media, Int. J. Numer. Methods Heat Fluid Flow (2019), https://doi.org/10.1108/HFF-022019-0129. [37] Y. Wang, S.Y. Sun, B. Yu, Acceleration of gas flow simulations in dualcontinuum porous media based on the mass-conservation POD method, Energies 10 (2017) 1380. [38] G. Jansen, B. Valley, S.A. Miller. THERMAID – a matlab package for thermohydraulic modeling and fracture stability analysis in fractured reservoirs, 2018. [39] D.W. Peaceman, Interpretation of well-block pressures in numerical reservoir simulation, Soc. Petrol. Eng. J. 18 (3) (1978) 183–194. [40] D.W. Peaceman, Interpretation of well-block pressures in numerical reservoir simulation with nonsquare gridblocks and anisotropic permeability, Soc. Petrol. Eng. J. 8 (3) (1983) 183–194. [41] K.A. Lie. An Introduction to Reservoir Simulation Using MATLAB: User Guide for the Matlab Reservoir Simulation Toolbox (MRST). SINTEF ICT, December, 2016 [42] M. Karimi-Fard, L.J. Durlofsky, K. Aziz, et al., An efficient discrete fracture model applicable for general purpose reservoir simulators, SPE Reservoir Simulation Symposium, Society of Petroleum Engineers, 2003. [43] S.B. Pluimers, Hierarchical Fracture Modeling Approach Master Thesis, Delft University of Technology, 2015. [44] L. Sirovich, Turbulence and the dynamics of coherent structures. I. Coherent structures, Q. Appl. Math. 45 (3) (1987) 561–571. [45] L. Sirovich, Turbulence and the dynamics of coherent structures. II. Symmetries and transformations, Q. Appl. Math. 45 (3) (1987) 573–582. [46] L. Sirovich, Turbulence and the dynamics of coherent structures. III. Dynamics and scaling, Q. Appl. Math. 45 (3) (1987) 583–590. [47] Y. Wang, B. Yu, S.Y. Sun, Fast prediction method for steady-state heat convection, Chem. Eng. Technol. 35 (4) (2012) 668–678. [48] Y.C. Zeng, J.M. Zhan, N.Y. Wu, et al., Numerical simulation of electricity generation potential from fractured granite reservoir through vertical wells at Yangbajing geothermal field, Energy 103 (1) (2016) 290–304. [49] T. Hayase, J.A.C. Humphrey, R. Greif, A consistently formulated quick scheme for fast and stable convergence using finite-volume iterative calculation procedures, J. Comput. Phys. 94 (1) (1991) 252.