A coupled dual continuum and discrete fracture model for subsurface heat recovery with Thermoporoelastic effects Текст научной статьи по специальности «Физика»

Математические заметки СВФУ Январь—март, 2019. Том 26, № 1

UDC 519.63

A COUPLED DUAL CONTINUUM AND DISCRETE FRACTURE MODEL FOR SUBSURFACE HEAT RECOVERY WITH THERMOPOROELASTIC EFFECTS

D. A. Ammosov, M. V. Vasilyeva, M. Babaei, and E. T. Chung

Abstract. We consider heat recovery from geothermal fractured resources with thermo-poroelastic effects. To this end, a hierarchical fracture representation is considered, where small-scale highly connected fractures are represented by the classical dual porosity model whereas large scale dense fractures are represented by the discrete fracture model. The mathematical model is described by a system of equations for mass and heat transfer for coupled dual continuum model as well as discrete fractures. Geomechanical deformations are written in the general form. For numerical solution of the resultant coupled system of equations including multicontinuum temperatures, pressures and deformations, we use the finite-element method. Numerical results are presented for two- and three-dimensional examples, showing applicability of the proposed method.

DOI: 10.25587/SVFU.2019.101.27250

Keywords: thermoporoelasticity, heat recovery, double porosity and double permeability, dual continuum, discrete fracture model, finite element method, mathematical modeling.

1. Introduction

Thermoporoelastic mathematical models are used for simulation of heat recovery from Enhanced Geothermal Systems (EGS) with applications on nuclear waste disposal and other coupled subsurface problems [1—6]. The mathematical model is described by coupled system of equations for mass transfer, heat transfer and geomechanical deformations of the porous media [7—11]. For simulation of the fractured porous media with dense distribution of the highly connected fractures, dual porosity models are used [12-15]. In dual porosity models, we suppose that fracture and porous matrix continua exist in each mesh cell.

For the cases when fracture distribution is sparse, we cannot suppose that fracture exists in each grid cell. For the explicit consideration of the fracture's geometry a mixed dimensional problem is used to describe a complex interaction between fracture and matrix. We can use embedded or discrete fracture models, that differ by a mesh construction for porous matrix and fracture. In the embedded

D. A. Ammosov's and M. V. Vasilyeva's works are supported by the mega-grant of the Russian Federation Government (No. 14.Y26.31.0013) and Russian Science Foundation (No. 17-71-20055).

© 2019 D. A. Ammosov, M. V. Vasilyeva, M. Babaei, and E. T. Chung

fracture model, mesh for fracture network and porous matrix is constructed separately, but for discrete fracture model, we construct conforming fracture mesh with porous matrix mesh [16-19].

In the numerical simulations, we should take into account all heterogeneities of the porous media, for example, in the geothermal reservoirs we can have different scales for fracture networks, one is for highly connected natural fracture networks and another for the large scale hydraulic fractures [20-23]. For describing the flow and heat transfer in the natural fracture networks, we can use dual continuum model, where small scale fractures and porous media matrix exist in each element of the grid mesh. For simulation of the large scale hydraulic fractures, we will use a discrete fracture model [19, 24-26]. Mechanic deformation is described using multicontinuum approach, where effective stress contains parts for each continuum for pressure and temperature [27-29].

The structure of the paper is as follows. In Section 2, we present mathematical model for the thermoporoelasticity problem in fractured porous media. Approximation of the model using a finite element method is presented in Section 3. In Section 4, we present numerical results for two and three-dimensional problems and finally, we present conclusions in Section 5.

2. Mathematical model

Here we present the thermoporoelasticity problem in fractured porous media. The mathematical model is described by the coupled system of equations for the deformations/displacements u, pressure p, and temperature T.

(a) To capture the small scale connected network of the fractures, we use a double porosity and double permeability (dual continuum) model. We introduce the two continua in the domain £ C Rd (d = 2, 3): the first is the small scale connected fracture network and the second is the porous matrix. Two continua are coupled by a transfer term. For each continuum, we define corresponding pressure and temperature, p1,T1 and p2,T2, respectively.

(b) To account large-scale fractures, we introduce domain 7 C Rd-i and define pressure pf and temperature Tf on 7 that connected by the mass and heat transfer terms.

For displacement u G £ C Rd, we use Biot's model with effective stress tensor that depends on pressure and temperature of each continuum. As such, we have the following coupled system of equations for pressure, temperature and effective displacements in three continua:

• Pressures (p1,p2,pf):

dpi dTi dV- u , ^ Sl-7^- - Vl-Qj- + ai~Ql--V(fciVpi) + 0-12 (pi -P2) + <71/(Pi - Pf) = 0,

dp2 0T2 dV- u , . S2 -Qj- - V-i-Qj- + a2 —Qt--V(fc2Vp2) + 021 (P2 -p 1) + 02/ (P2 -P}) = 0, UJ

dpf dTf

sf— -Vf-gf - Vifc/Vp/) + 0/i(p/ - Pi) +072(p/ -p2) = 0,

where si = 1/Mi; Mi is the Biot modulus (i = 1, 2, f ). • Temperatures (T1,T2,Tf ):

^ , nwyitrr a , « rr dV ■ u dpi

Cl~gf + Ci V(Tl 'QV ^ lfi~dt--V '°~dt

-V(AiVTi)+ r 12(Ti - T2)+ rif (Ti - Tf) = 0,

^ , nwv7lrr x , « ^ dV • u dp2 + 2 V(T2 • Ç2) + /32^2,0—^--^T^fi —

-V(A2VT2) + r2i(T2 - Ti) + r2f (T2 - Tf ) = 0,

(2)

dTf

dpf

dt

-V(Af VTf ) + rf i(Tf - Ti) + rf2(Tf - T2) = 0, where qi = -kiVpi is the Darcy's velocity, ki = Kj/vf is the mobility of the i-th continuum , Kj is the permeability of the i-th continuum, Vf is the fluid viscosity (i = 1, 2, f ).

• Displacements u:

-Va + "jVPi + X) &VTi = 0, (3)

i=i,2,f i=i,2,f where a and e are the stress and strain tensors

o- = Atr(e)/ + 2yue, e = -(Vit + (Vw)T),

and A, y are the Lame coefficients.

Here Aj is the thermal conductivity of the i-th continuum, ai is the Biot-Willis constant, pi = ^i(A + 2/3«), ^ is thermal expansion of the solid phase, Ci is the equivalent volumetric heat capacity, n is the thermal expansion coefficient, aij and rij are the mass and heat transfer coefficients [30].

We supplement system of equations (1), (2) and (3) with the following initial conditions

Pi = Pi,0, Ti = Tj,o, u = uo (4)

and the boundary conditions:

• for pressure:

dpi

~kl~Én =xGrj' ier¿.

=0, x G dil,

dn

dpf

—kf—— = 0, x G r¿v, pf = g, x G Ir>, dn

(5)

where TJ U TL = and TD U TN = dy; • for temperature:

dTi

—Ai—— = 0, x G Tj, Ti = h, i £ dn

dT

-A2^ = 0, a; G dn

dTf

—Af—— = 0, x G Tat, Tf = h, x G l£>; dn

• for displacement:

a ■ n = 0, x £ 11, u = 0, x £ r2, (7)

where Ti U T2 = dQ.

The above triple continuum model can be written in the following general form of the multicontinuum model

dpi dTi d V ■ u

Sl~dt ~ rll~dt + ai~dt--Vi^vp^ + 2_^o-ij{Pi ~Pj) = 0,

j

n dTl i nwyj/rr \ , a rr dV ' u rr dpi

dt dt dt (g)

- V(AiVTi) + ^ rij(Ti - Tj) = 0,

j

-Va + Y, aVPi + ^2 frVTi = 0, ii

where qi = -kiVpi, pi and Ti are the pressure and temperature of the i-th continuum, where i = 1, 2,.... We note that, in this model we ignore gravity forces and suppose diffusive transfer term for the heat equations.

3. Finite-element approximation

For numerical solution of the problem, we use a finite element method with discrete fracture model (DFM). We construct an unstructured grid, where fractures are located on the interfaces between the cells. We use the continuous Galerkin finite element method with linear basis functions, and by using the superposition principle, we eliminate equation for the fracture network.

For finite element approximation, we define space V = H(Q), Vg = {v £ H(Q) : vk = g}, V7,g = {v £ H(y) : vliD = g}, W = {w £ [H(Q)]d : w|r2 = 0} and u £ W, pi £ Vg, Ti £ Vh, P2 £ V, T2 £ V, pf £ V7,g, Tf £ V7,h. We have following variational formulation of the problem: find (p1,p2,pf, T1,T2,Tf ,u) £ Vg x V x VYg x Vh x V x V7ih x W, such that

f pi - pi , f Ti - Ti ,

/ si-t>i ax — / 771-t>i ax

J T J T

n n

r V ■ u - V ■ u , /•,__,

+ / ai-t>i ax + / kiVpi ■ vvi ax

n n

+ j ai2(pi - p2)vi dx + / aif (pi - pf )vi dx = 0, (9)

P2 - p2 , f T2 - T2 S 2-V2dx — / 772-1>2 UX

n

r V-u-V-u , /•,__,

+ / «2-V2dx + / K2 Vp2 ' Vf2 UX

n n

— J ^12 (Pi - P2)v2 dx + J 0*2/ (p2 - P/ )V2 dx = 0, (10)

P / - P / f T/ - Í1/ f

Sf—--Vfdx— / í//—--Vf dx + / kf'Vpf -S7vf dx

-j ^ (pi - p/)v/ dx-J ^(p2 - p/)v/ dx = 0, (11)

J CiTl t Tl zi ¿E + J C7V(Ti-gi)zi<ic

nn

f V-u-V-u f pi - pi

+ / pi J 1,0-zi ax — / rjilifi-zi dx

n

+ J A1VT1 ■ Vzi dx + J ri2(Ti - T2)zi dx + J ri/(Ti - T/)zi dx = 0, (12)

n n Y

y C2T2 ~T2g2dx + y C™V{T2-q2)z2dx n n

[ V • U - V • U f P2 p2 ,

+ / P2-12,0-z2dx- / 7722,0-z2ax

n

+ y A2VT2 ■ Vz2 dx - J ri2 (Ti - T2)z2 dx + y r2/ (T2 - T/)z2 dx = 0, (13) n n Y

y cf Tf ~ T/ + y CfV(Tf<lf)zfdx- J VfTf/f ~pf Z2 dx

Y Y n

+ / A/VT/ ■ Vz/ - / ri2m - T/)„ dx - / ,,/ffi - T,)„ dx = 0 <14>

Y Y Y

y o-(u) : e(w) dx + j aiVpi ■ wdx + j a2Vp2 ■ wdx n n n

+ JßiVTi ■ wdx + y ß2VT2 ■ wdx = 0, (15)

where (vi,v2,vf ,zi,z2, zf ,w) G V0 x V x VY$ x V0 x V x VY,0 x W. We note that, we used an implicit scheme for approximation by time with time step t, r is the solution from previous time step.

Let ¿Ph be a mesh for the domain £2, £\ be the faces of the grid ¿Ph, and SY will be the faces of large-scale fractures, where S'Y C £\. Using the principle of superposition, we obtain the following discrete system in the matrix form for y = (pi,P2,Ti,T2,u)T

where

A =

For matrices, we have • Si = {si,lj + Lisf,ij}, I =1, 2,

My-y + Ay = 0,

t

S1 0 Ni 0 Gi \

0 S2 0 N2 G2

M = N 0 Hi 0 Bi ,

0 N2 0 H2 B2

00 0 0 0

+ Q -Q 0 0 0

-Q D2 + Q 0 0 0

0 0 Ai + R -R 0

0 0 -R A2 + R 0

GT G2T BT B2T kJ

sMj = / ^ dx, s2i = j S2^3 dX, f = / sf dX,

n n Y

• Hi = {ci,ij + licf,ij}, I =1, 2,

^ = / dX' C2i = / dX, Cfi = / Cf *f>f dX' n n Y

• Ni = TifiNi, Ni = {nij + linfj}, I =1, 2,

m« = -J niMj dx, n2,j = -J V2Mj dx, f = -J Vf dx,

n n Y

• Di = {di,ij + iidf,ij}, I = 1, 2,

di,j = / kM-V^j dx d2,j = j kM-V^j dx dfi =j kf dx

n n Y

• Ai = {ai,ij + iiaf,ij}, I = 1, 2,

a1,ij = j A1V^i • V^j dx + j CW• qi)^j dx,

a2,ij = J A2V^i • V^j dx + J CW• q_2)^j dx

af,ij = A/Vf • Vf dx + CJV(f • qf dx,

• K = {kj} with kij = J a(^) : e(®j) dx,

Q

• Gi = {gij} with gij = f ai V • u($i dx, g2j = f V • u($i dx,

Q Q

l = 1, 2,

• Bi = {bi,ij} with biij = / pi V • u($idx, b2ij = f p2V • u($idx,

QQ

l = 1, 2,

• Q = Uij } and R = {rij } with qj = f a 12 fai faj dx, rj = f ri2 ^i ^j dx.

QQ We use linear basis functions for domain Q and lower dimensional basis for domain 7. We use artificial diffusion for stabilization of the temperature equation, Ai = Xl + Xad (ql ) and take velocity from the previous time step pressure solution. Here, we eliminate equations for the fractures by superposition principle and by supposing equality of the pressure and temperature at the fracture nodes.

4. Numerical results

In this section, we present numerical results for the thermoporoelasticity problem defined using presented finite element approximation. We consider two test problems:

• Two-dimensional problem in domain Q = [0,10] x [0,10] m2,

• Three-dimensional problem in domain Q = [0,10] x [0,10] x [0, 5] m3.

Numerical implementation is based on the open source library FEniCS [31]. In order to resolve fracture networks by the mesh, we construct an unstructured grids using Gmsh software.

Fig. 1. Computational grids for two- and three-dimensional problems. Left: two-dimensional domain Q. Middle: three-dimensional domain Q. Right: fracture networks 7

We set following parameters for the computational grids:

• Two-dimensional problem: 3944 vertices and 7685 cells,

• Three-dimensional problem: 23664 vertices and 115275 cells.

1 000e+07 2.000e+07 1.000e+07 2.000e+07 1.000e+07 2.000e+07

Fig. 2. Pressure distribution (Pa) for different times tn, n = 5, 25 and 100

(from left to right) for two-dimensional test case. The first row: the first continuum. The second row: the second continuum.

Fig. 3. Temperature distribution (°C) for different times tn, n = 5, 25 and 100 (from left to right) for two-dimensional test case. The first row: the first continuum. The second row: the second continuum.

-0.0138 -0.0069 0.0000 -0.0043 -0.0024 0.0005

Fig. 4. Displacements distribution (m) at final time for two-dimensional test case, ux and uy (from left to right).

1.000e+07 2.000e+07 1.000e+07 2.000e+07 1.000e+07 2.000e+07

Fig. 5. Pressure distribution (Pa) for different times tn, n = 5, 25 and 100 (from left to right) for three-dimensional test case.

The first row: the first continuum. The second row: the second continuum

We depict grids in Fig.1, where we represent unstructured grid with triangular cells for two-dimensional domain Q on the left picture, an unstructured grid with tetrahedral cells for the first and the second continuum in the domain Q and fracture domain y are represented in the middle and right pictures for three-dimensional model problem.

We use following parameters in numerical simulations:

• Fluid flow:

S1 = S2 = 10-7, Sf = 0 [Pa-1 ], ni = n2 = 3 • 10-4, f = 0[°C-1 ],

ki = 10-12, k2 = 10-13, kf = 10-9 [m2 • Pa-1 • s-1 ],

• Heat transfer:

1.0006407 2.000e+07 1.000e+07 2,000e+07 l.OOOe+07 2 000e+07

Fig. 6. Temperature distribution (°C) for different times tn, n = 5, 25 and 100 (from left to right) for three-dimensional test case. The first row: the first continuum. The second row: the second continuum

1.0006407 1.004e+07 1.000e+07 1.019e+07 l.OOOe+07 1 073e+07

Fig. 7. Displacements distribution (m) at final time for three-dimensional test case, ux, Uy and uz (from left to right)

C1 = C2 = 107, Cf = 0, [J • m-3 • °C-1] A1 = 100, A2 = 10, Af = 105 [W • m-1 • °C-1],

• Geomechanical properties: a1 = a2 = 0.1, af = 0,

£1 = £2 = 3 • 10-5(A + 2/3y), £f = 0 [°C-1 • Pa], A = m = 109 [Pa].

Here, we suppose that (1) the first continuum represents the dense small scale highly connected fracture networks, (2) the second continuum represents the porous matrix, and (3) the third continuum represents the large scale fracture network. For

transfer terms, we set <r12 = 10-13 and r12 = 10 for interaction between the first and the second continuum, and <r2f = r2f = 0, that means that the third continuum interacts only with the first continuum.

We simulate for tmax = 86400 [s] with 100 time steps. We have following initial and boundary conditions:

• Initial conditions: T0 = 200 [°C] and p0 = 107 [Pa],

• Boundary conditions: h = 20 [°C] and g = 2 • 107 [Pa] on the left boundary for the first and the third continuum. On other boundaries, we set zero Neumann boundary conditions.

Results of the numerical simulations are presented in Fig. 2, 3, and 4 for two-dimensional test case and in Fig. 5, 6, and 7 for three-dimensional test. In Fig. 2 and 5, we depict pressures for the first and the second continua for different time steps, ¿5, t25 and t1oo. Temperature distributions for the first and the second continua at times of ¿5, t25 and t1oo are presented in Fig. 3 and 6. Displacement distribution at the final time are presented in Fig. 4 and 7 for two- and three-dimensional test problems, respectively. We observe similar behavior for the two and three-dimensional test cases. We observe how pressure increase from left to right, when we inject cold water. Pressure differences lead to the temperature propagation from the left boundary. For the first continuum (small scale fractures), we can see large pressure and temperature drop due to injection. For the second continuum (porous matrix), we observe smaller gradients, because we set <r2f = r2f = 0 in our hierarchical model. This means that we have heat and mass flow from large scale fractures (the third continuum) to small scale fractures (the first continuum) and after that the first continuum interacts with the second continuum. In Fig. 4 and 7, we observe deformations that are caused due to pressure and temperature gradients (the results have been warped by displacements vector multiplied by 120). We note that, the model that we consider in this work contains several simplifying assumptions. We will consider full model in the future works with some multiscale model reduction techniques.

5. Conclusion

We considered the thermoporoelasticity problem in fractured media and presented a general multicontinuum model for mass and heat transfer for a coupled dual continuum model and discrete fractures to determine geomechanics deformations. We constructed a finite-element approximation for the coupled system of equations for pressure, temperature and displacements. Numerical results are presented for the two and three-dimensional formulations.

REFERENCES

1. Karvounis D. C., Simulations of Enhanced Geothermal Systems with an Adaptive Hierarchical

Fracture Representation, Diss. Doct. Sci., ETH Zurich, Zurich (2013).

2. Sandve T. H., Berre I., Keilegavlen E., and Nordbotten J. M., "Multiscale simulation of flow

and heat transport in fractured geothermal reservoirs: inexact solvers and improved transport

upscaling," in: 38th Workshop Geothermal Reservoir Engineering 2013, Stanford Geothermal Program, 1, 575-583 (2013).

3. Akkutlu I. Y., Efendiev Y., and Vasilyeva M., "Multiscale model reduction for shale gas transport in fractured media," Comput. Geosci., 20, No. 5, 953-973 (2016).

4. Akkutlu I. Y., Efendiev Y., Vasilyeva M., and Wang Y., "Multiscale model reduction for shale gas transport in a coupled discrete fracture and dual-continuum porous media," J. Nat. Gas Sci. Eng., 48, 65-76 (2017).

5. Praditia T., Helmig R., and Hajibeygi H., "Multiscale formulation for coupled flow-heat equations arising from single-phase flow in fractured geothermal reservoirs," Comput. Geosci., 22, No. 5, 1305-1322 (2018).

6. Salimzadeh S., Nick H. M., and Zimmerman R. W., "Thermoporoelastic effects during heat extraction from low-permeability reservoirs," Energy, 142, 546-558 (2018).

7. Wang X. and Ghassemi A., "A 3D thermal-poroelastic model for geothermal reservoir stimulation," in: 37th Workshop Geothermal Reservoir Engineering 2012, Stanford Geothermal Program, 1, 936-946 (2012).

8. Karvounis D. C. and Jenny P., "Adaptive hierarchical fracture model for enhanced geothermal systems," Multiscale Modeling & Simulation, 14, No. 1, 207-231 (2016).

9. Hao Y., Fu P., and Carrigan C. R., "Application of a dual-continuum model for simulation of fluid flow and heat transfer in fractured geothermal reservoirs," in: 38th Workshop on Geothermal Reservoir Engineering 2013, Stanford Geothermal Program, 1, 462-469 (2013).

10. Nissen A., Keilegavlen E., Sandve T. H., Berre I., and Nordbotten J. M., "Heterogeneity preserving upscaling for heat transport in fractured geothermal reservoirs," Comput. Geosci., 22, No. 2, 451-467 (2018).

11. Salimzadeh S., Paluszny A., Nick H. M., and Zimmerman R. W., "A three-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal systems," Geothermics, 71, 212-224 (2018).

12. Barenblatt G. I., Zheltov Iu. P., and Kochina I. N., "Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]," J. Appl. Math. Mech., 24, No. 5, 1286-1303 (1960).

13. Arbogast T., Douglas J., Jr., and Hornung U., "Derivation of the double porosity model of single phase flow via homogenization theory," SIAM J. Math. Anal., 21, No. 4, 823-836 (1990).

14. Kazemi H., Merrill L. S., Jr., Porterfield K. L., and Zeman P. R., "Numerical simulation of water-oil flow in naturally fractured reservoirs," Soc. Petrol. Eng. J., 16, No. 6, 317-326 (1976).

15. Warren J. E. and Root P. J., "The behavior of naturally fractured reservoirs," Soc. Petrol. Eng. J., 3, No. 3, 245-255 (1963).

16. Vasilyeva M., Chung E. T., Leung W. T., and Alekseev V., "Nonlocal multicontinuum (NLMC) upscaling of mixed dimensional coupled flow problem for embedded and discrete fracture models," 2018. arXiv preprint arXiv:1805.09407

17. Karimi-Fard M., Durlofsky L. J., and Aziz K., "An efficient discrete fracture model applicable for general purpose reservoir simulators," SPE Reservoir Simulation Symp., Soc. Petrol. Eng., 2003.

18. Lee S. H., Jensen C. L., and Lough M. F., "Efficient finite-difference model for flow in a reservoir with multiple length-scale fractures," Soc. Petrol. Eng. J., 5, No. 3, 268-275 (2000).

19. Vasilyeva M., Chung E. T., Cheung S. W., Wang Y., and Prokopiev G., "Nonlocal multicon-tinua upscaling for multicontinua flow problems in fractured porous media," arXiv preprint arXiv: 1807.05656 (2018).

20. Martin V., Jaffre J., and Roberts J. E., "Modeling fractures and barriers as interfaces for flow in porous media," SIAM J. Sci. Comput., 26, No. 5, 1667-1691 (2005).

21. D'Angelo C. and Quarteroni A., "On the coupling of 1D and 3D diffusion-reaction equations: application to tissue perfusion problems," Math. Models Meth. Appl. Sci., 18, No. 8, 14811504 (2008).

22. Formaggia L., Fumagalli A., Scotti A., and Ruffo P., "A reduced model for Darcy's problem in networks of fractures," ESAIM: Math. Model. Numer. Anal., 48, No. 4, 1089-1116 (2014).

23. D'Angelo C. and Scotti A., "A mixed finite element method for Darcy flow in fractured porous media with non-matching grids," ESAIM: Math. Model. Numer. Anal., 46, No. 2, 465-489 (2012).

24. Gong B., Karimi-Fard M., and Durlofsky L. J., "Upscaling discrete fracture characterizations to dual-porosity, dual-permeability models for efficient simulation of flow with strong gravitational effects," Soc. Petrol. Eng. J., 13, No. 1, 58-67 (2008).

25. Chung E. T., Efendiev Y., Leung T., and Vasilyeva M., "Coupling of multiscale and multi-continuum approaches," GEM-Int. J. Geomath., 8, No. 1, 9-41 (2017).

26. Vasilyeva M., Babaei M., Chung E. T., and Spiridonov D., "Multiscale modelling of heat and mass transfer in fractured media for enhanced geothermal systems applications," Appl. Math. Model., 67, 159-178 (2018).

27. Vasilyeva M. and Tyrylgin A., "Machine learning for accelerating effective property prediction for poroelasticity problem in stochastic media," arXiv preprint arXiv:1810.01586 (2018).

28. Vasilyeva M., Chung E. T., Efendiev Y., and Kim J., "Constrained energy minimization based upscaling for coupled flow and mechanics," arXiv preprint arXiv:1805.09382 (2018).

29. Akkutlu I. Y., Efendiev Y., Vasilyeva M., and Wang Y., "Multiscale model reduction for shale gas transport in poroelastic fractured media," J. Comput. Phys., 353, 356-376 (2018).

30. Coussy O., Poromechanics, John Wiley & Sons, Chichester (2004).

31. Logg A., Mardal K.-A., and Wells G. (ed.), Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Springer Science & Business Media, New York (2012).

Submitted November 29, 2018 Revised December 28, 2018 Accepted March 1, 2019

Dmitry A. Ammosov Multiscale model reduction laboratory, North-Eastern Federal University, Yakutsk 677980, Russia dmitryammosov@gmail.com

Maria V. Vasilyeva

Department of Computational Technologies,

North-Eastern Federal University,

Yakutsk 677980, Russia;

Institute for Scientific Computation,

Texas A&M University,

College Station, TX 77843-3368

vasilyevadotmdotv@gmail.com

Masoud Babaei

The University of Manchester,

School of Chemical Engineering and Analytical Science, Manchester, M13 9PL, UK masoud.babaei@manchester.ac.uk

Eric T. Chung

Department of Mathematics,

The Chinese University of Hong Kong (CUHK),

Hong Kong SAR

tschung@math.cuhk.edu.hk