A sequential implicit discrete fracture model for three-dimensional coupled flow-geomechanics problems in naturally fractured porous media

A sequential implicit discrete fracture model for three-dimensional coupled flow-geomechanics problems in naturally fractured porous media

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx Contents lists available at ScienceDirect Journal of Petroleum Science and Engineering...

2MB Sizes 1 Downloads 95 Views

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol

A sequential implicit discrete fracture model for three-dimensional coupled flow-geomechanics problems in naturally fractured porous media ⁎

Mostafa Moradia, Amir Shamlooa, , Alireza Daneh Dezfulib a b

Department of Mechanical Engineering, Sharif University of Technology, Azadi Ave., Tehran, Iran Department of Mechanical Engineering, Faculty of Engineering, Shahid Chamran University, Ahvaz, Iran

A R T I C L E I N F O

A BS T RAC T

Keywords: Fractured porous media Discrete-fracture model Coupled flow-geomechanics Finite Element method

A sequential implicit numerical method based on discrete-fracture model and the Galerkin Finite Element method, for time-dependent coupled fluid flow and geomechanics problems in fractured subsurface formations is presented. Discrete-fracture model has been used to explicitly represent the fracture network inside porous media. The Galerkin Finite Element method with adaptive unstructured gridding is implemented to numerically solve the spatially three-dimensional and time-dependent problem. The presented method is validated with previously obtained solutions. Two problems are numerically solved by applying the presented methodology in a three-dimensional fractured petroleum reservoir under different production and injection conditions. The modeling results demonstrate that changes in fracture aperture and permeability due to geomechanical stress, have a significant impact on well performance and production rate in fractured reservoirs.

1. Introduction The stress dependent fluid flow is one of the challenging problems in many fields of geosciences. Interactions between solid media and the fluid which is stored in it, lead to an inherently time-dependent coupled fluid flow and geomechanics problem. These interactions are due to changes in fluid pore pressure (via injection or removal of the fluid), which in turn cause the solid media to deform. As a result of this deformation, the characteristics of the solid media such as permeability (i.e. the ability of the media to conduct fluid) change. The changes in solid media's characteristics alters the way that the fluid flows within it. This chain of events inside the porous media demonstrates the nature of the coupled fluid flow and geomechanics problem (Rutqvist and Stephansson, 2003). Disregarding the effects of geomechanical parameters on fluid flow in petroleum reservoirs can lead to significant economic and environmental damages. For instance, excessive water flooding for oil recovery enhancement, without considering the geomechanical factors and how they activate or deactivate fluid pathways has caused massive oil leaks in the past (Gu et al., 2014). Therefore it is important to pay careful attention to geomechanical effects in fractured formations. Many previous studies exist in the literature, concerning the stresssensitivity effects in conventional porous media (Pedrosa, 1986; Kilmer et al., 1987; Zhang and Ambastha, 1994; Chin et al., 2000; Wang et al., 2010; Qanbari and Clarkson, 2012, 2014; Shaoul et al., 2015).



However, coupled fluid flow and geomechanics study in fractured formations is a relatively new field of study in reservoir simulation. Different approaches have been taken by several researchers to present a method for solving the stress dependent fluid flow problem in fractured porous media. In one-way coupled method, the changes in fluid properties affect the geomechanics, but geomechanical information is not shared with the fluid flow. On the other hand, two-way coupled methods allow fluid flow and geomechanics to interact with each other. Numerical simulation of two-way coupled problems can be carried out using two distinct approaches: one is the fully implicit scheme, in which the flow and geomechanics equations are solved simultaneously, and the second approach is called the sequential method. Fully implicit methods offer a more accurate solution, but have large computational costs. However, sequential methods are more flexible from a computational perspective, but might suffer from convergence and stability limitations (Castelletto et al., 2015; Gutierrez et al., 2001; Kim et al., 2012; Jalali and Dusseault, 2012; Settari and Mourits, 1998). Another approach which is called loose coupling, is a middle ground between one-way and two-way coupling method. In loose coupling, two sets of equations, one for fluid flow and the other for geomechanics are solved independently, but information is shared between the two at specified time intervals (Minkoff et al., 2003). There are different numerical methods to solve the fluid flow and the geomechanics problem in porous media. For the fluid flow problem,

Corresponding author. E-mail address: [email protected] (A. Shamloo).

http://dx.doi.org/10.1016/j.petrol.2016.12.027 Received 13 September 2016; Received in revised form 20 November 2016; Accepted 20 December 2016 0920-4105/ © 2016 Elsevier B.V. All rights reserved.

Please cite this article as: Moradi, M., Journal of Petroleum Science and Engineering (2016), http://dx.doi.org/10.1016/j.petrol.2016.12.027

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

finite element based methods have been widely used to simulate the coupled problem (Karimi-Fard and Firoozabadi, 2003; Monteagudo and Firoozabadi, 2004; Martin et al., 2005). Many other works based on the finite volume method have also been published (Karimi-Fard, 2004). For the solid mechanics part, finite element based methods have also been used extensively (Zienkiewicz and Taylor, 2005). Finite element based methods and hybrid finite element-finite volume methods have been employed in solving the coupled fluid flow and geomechanics problems (Gu et al., 2014; Garipov et al., 2016). Integral finite difference method is another approach that has been used in such problems (Hu et al., 2013). Another relatively new scheme, is the extended finite element method (XFEM) which is suitable for problems with local features that cannot be fixed by mesh refinement (Khoei et al., 2006; Lamb et al., 2013; Shamloo et al., 2005). Many researchers have investigated the coupled problem in fractured media. Some of these works have utilized the dual-porosity model (Chen and Tuefel, 2000; Segura et al., 2016). Discrete fracture model has been more widely used recently in coupled fluid flow and geomechanics problems (Garipov et al., 2016; McClure et al., 2016; Norbeck et al., 2014; Philip et al., 2005). Instead of solving the coupled system of equations, Moinfar et al. (2013) used in-situ stress conditions to estimate the stress field within the reservoir, and utilized an experimental relation to compute the fracture aperture deformation due to the stress field (Moinfar et al., 2013). Rutqvist et al. (2014) implemented a linked multicontinuum and crack tensor approach to simulate coupled fluid flow, geomechanics and solute transport in fractured rock (Rutqvist et al., 2014). Garipov et al. (2016) presented a fully implicit coupled fluid flow and geomechanics formulation for fractured formations. In their work, contact problem has been included to study the effects of geomechanical stress on fracture aperture and hydrocarbon production rate (Garipov et al., 2016). In the present work, we aim to implement the discrete fracture model along with the Galerkin finite element model to present a sequential implicit time-dependent numerical method to solve the coupled fluid flow and geomechanics problem in subsurface naturally fractured formations. The presented method is sequential implicit, which means that first, the fluid flow equation is solved implicitly. After obtaining the unknown fluid flow variable, i.e. pore pressure, geomechanical equations are explicitly solved based on the values acquired for pore pressure. We use the discrete fracture model (Karimi-Fard and Firoozabadi, 2003; Monteagudo and Firoozabadi, 2004) to represent the fracture system inside the porous media. In this work, an experimental relation is used to correlate the fracture aperture deformation with the calculated stress field. Firstly, the theory and governing equations for fluid flow and geomechanical stress, fracture modeling and fracture aperture relation with stress field are presented. In the third section, the numerical method and computational grid structure are explained. The method is then validated with the results obtained from similar research in section four. After model validation, some examples of the application of the presented method are offered.

Fig. 1. Shear dilation.

• •

In-situ stress field is constant at all times. Gravity effects are neglected.

2.2. Governing equations The mechanical behavior of Porous media is described by the poroelasticity theory. According to this theory, the total stress is the sum of stress caused by fluid pore pressure and the stress inside the rock skeleton (Biot, 1941; Chen et al., 1995).1

σij =σij′+αPδij

⎧ 1,& i = j δij=⎨ ⎩ 0, & i ≠ j

Where σ is the total stress tensor, σ′ is the biot effective stress, P is fluid pore pressure and α is the biot coefficient. Fluid flow equation is governed by the law of conservation of mass:

α

∂ϵ ∂P +S =−∇.v ∂t ∂t

(2)

Where ϵ is strain, v is fluid velocity and S is storativity which is equal to:

S =ϕCf +(α −ϕ) Cr

(3)

Where, ϕ is the porosity of porous media, Cf is fluid compressibility and Cr is rock skeleton compressibility. In Eq. (2), there are two unknown variables, P and v . However according to Darcy's law, fluid velocity is related to pore pressure:

k v=− ∇P μ

(4)

Where k is the permeability of porous media and μ is the viscosity of fluid. Combining, Eqs. (3) and (4), results in Eq. (5).

α

∂ϵ ∂P +S ∂t ∂t

k =∇.( ∇P ) μ

(5)

The stress which acts on a representative volume of the porous media is the total stress. The conservation of momentum in rock skeleton can be written as:

∇.σ −f =0 2. Mathematical formulation

(6)

Where f is the vector of external forces. Assuming that deformations are small and porous media is linearly elastic, effective stress can be written in terms of strain:

In this section, the underlying theory and equations governing the coupled fluid flow and geomechanics problem are presented. The assumptions made in the method are also presented.

⎛ 2 ⎞ σij′=−⎜K − G⎟ ϵ v −2G ϵij i=jσij′=−2G ϵij i=j ⎝ 3 ⎠

2.1. Assumptions

• • •

(1)

(7)

In Eq. (7), K is the bulk modulus and G is the shear modulus of rock skeleton, ϵ v is the volumetric strain and is equal to:

The fluid is single-phase and slightly compressible and has constant properties. The solid media is linearly elastic and its properties also remain constant. The fractures are two dimensional rectangular planes that have variable apertures that are functions of the applied stress field.

ϵ v=ϵxx +ϵ yy +ϵzz

1

2

Tensile stress is positive.

(8)

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

In Eq. (13), Kn is the normal fracture stiffness. 2- Shear stress, τf , acting along the fracture makes the fracture surfaces to slide relative to each other and the normal component of this displacement is the change in fracture aperture. This mechanism is called “shear dilation” and it is shown in Fig. 1. Shear dilation does not always occur when a shear stress is applied to the fracture. Shear dilation starts when the shear stress is equal to the shear strength, σc . In other words, when τf < σc , shear stress causes tangential displacement to occur between the fracture surfaces, however normal displacement does not happen until τf > σc (Barton et al., 1983). The relation between shear stress and tangential displacement is:

Strain tensor is related to rock skeleton displacement vector by Eq. (9).

ϵ=1/2(∇u+∇T u)

(9)

Where u is the rock skeleton displacement vector. Rewriting Eq. (6) in terms of the displacement vector and then combining it with Eq. (1) lead to the equation shown below:

⎛ 1 ⎞ ⎜K + G⎟ ∇ϵkk +G∇2 u−α∇P=f ⎝ 3 ⎠

(10)

So far, there are four independent variables (i.e. three components of displacement vector and pore pressure) and four equations (Eqs. (5) and (10)). Note that ϵkk is a function of u, and therefore is not an independent unknown variable.

dt =

τf (14)

Ks

2.3. Discrete fracture model

Nevertheless, the tangential displacement that is associated with shear dilation, is related to shear stress as follows:

Flow modeling in fractured formations is substantially more complicated than conventional porous media. The main difficulty associated with modeling fractured porous media is that fractures have a strong effect on the direction of flow and the way fluid is transported within the porous media. The difference between the size scale of fracture width and porous media renders simple numerical methods costly and difficult (Jaffre et al., 2011). Therefore it is necessary to implement a special method to overcome this difficulty. There are different approaches for modeling the system of fractures inside the porous media. Continuum models rely on averaging and homogenization methods, and therefore cannot model the geometry of the fracture system accurately. But on the hand, the discrete-fracture model represents each fracture individually. This allows us to include and model the system of fractures inside the porous media more accurately. In discrete-fracture model, the domain is subdivided into two domains: fracture system and porous media (matrix) and it is assumed that fractures are two-dimensional entities within a three-dimensional space. The system of Eqs. (10) and (5) are solved in these two subdomains (Monteagudo and Firoozabadi, 2004).

⎧ τ f − σc ,& τ > σ ⎪ f c dt , sd =⎨ Ks ⎪ ⎩ 0, & τf < σc

Ω = Ωm +df Ωf

dn, sd =dt , sd tan(Ψ)

(16)

In Eq. (16), Ψ is the shear dilation angle. The above mentioned relations are valid if the fracture is linearly elastic with respect to normal and tangential stresses applied to it. However, if stress levels exceed a certain value, plastic deformation occurs and different relation should be used to correlate displacement with the stress field. Using the Mohr-Coulomb criterion, fracture failure point can be determined. According to the Mohr-Coulomb criterion for rock media, maximum shear stress that can be applied to fracture, without failure is as follows (Olsson and Barton, 2001).

τf , max =σf tan(ϕ)

(17)

Where τf , max is the maximum shear stress before failure, and ϕ is the fracture internal friction angle. Assuming perfect plastic fracture behavior after failure, the following relations can be written for fracture aperture.

(11)

In Eq. (11) Ω represents the whole domain, m and f denote the matrix and fracture system, respectively, and d f is the fracture aperture. Since the fracture system,Ω f is two-dimensional and Ω m is a three-dimensional domain, in order to resolve the discrepancy between the dimensions, fracture system domain is multiplied by fracture aperture, d f . 2.4. Fracture aperture and mechanical stress As mentioned earlier, fracture system properties are functions of the reservoir mechanical stress field. Two distinct mechanisms cause the fracture aperture to change: 1- Variation of normal stress applied to fracture, directly alters the fracture aperture. The change in fracture aperture, Δdn is related to normal stress applied to the fracture, σf according to the well-known Barton and Bandis equation (Barton et al., 1985):

Δdn=

Dn 1+

Kn, i Dn σf

(12)

Where Δdn is the normal fracture closure, Dn is the maximum fracture aperture closure, and Kn, i is the initial normal fracture stiffness. Normal fracture stiffness is a measure of fracture's sensitivity to normal stress and can be calculated by Eq. (13).

⎡ ⎤−2 σf ⎥ Kn=Kn, i ⎢1− ⎣ Kn, i Dn + σf ⎦

(15)

Where, dt , sd is the tangential displacement associated with shear dilation, and Ks is the shear stiffness of fracture. The normal component of displacement caused by shear dilation is (Barton et al., 1985):

(13)

Fig. 2. Solution procedure in each time step.

3

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

d f =d f ,0−Δdn+min(dn, sd , Ucs ) dn, sd ⎧ 0 τf ≤σf ⎪ dt , sd tan(Ψ) ⎪ σ ≤ τ =⎨ ⎛ c f < τf , max ⎪ ⎜ τ f , max − σc ⎞⎟ tan(Ψ)+ σf − σf , max τf =τf , max ⎪ ⎝ Ks ⎠ Kn tan(Ψ) ⎩

(18)

Where σf , max is the normal stress corresponding to maximum shear stress, and Ucs is the maximum shear displacement. After fracture aperture is calculated, fracture permeability can also be calculated by Eq. (19).

kf =

d f2 (19)

12

Fig. 3. 2D validation problem.

3. Numerical method and solution strategy

Table 1 input data.

200m × 400m

In this section, the numerical method used to solve the system of Eqs. (5) and (10) is presented. Furthermore, solution strategy and procedure are also explained.

5 × 10−14m2 0.15

3.1. Galerkin Finite Element method

Rock properties Domain dimensions Matrix permeability, km Matrix Porosity, ϕ Compressibility, cr

6.29 × 10−10Pa−1 0.83

Biot’s coefficient, α Young’s modulus, E Poisson’s ratio, ν Internal friction angle, ϕ Shear dilation angle, Ψ Shear strength, σc

5.8 × 109Pa 0.3 24.9° 5

Initial normal stiffness, Kn, i

1.2 × 1011Pa /m

Shear stiffness, Ks

1.2 × 1011Pa /m

No-load fracture aperture, df , 0

2.5 × 10−4m

Maximum fracture aperture closure, Dn

2.4 × 10−4m

Fluid properties Viscosity, μ

1.8 × 10−4Pa. s

Compressibility, cf

5.39 × 10−9Pa−1

Initial and Boundary conditions Initial pore pressure, Pi

10 × 106Pa

Injecting well pressure, Pinj

20 × 106Pa

Producing well pressure, Pprod

10 × 106Pa

In-situ stress, σis

20 × 106Pa

In this work, we use the Galerkin variational method with Finite Element discretization to solve the system of Eqs. (5) and (10). Numerical methods based on the Finite Element scheme are numerically more demanding compared to other numerical methods, such as the Finite Difference method, but since they can be implemented in conjunction with unstructured gridding, Finite Element methods are suitable for three-dimensional problems with complex geometry. Another advantage of the Finite Element based methods is that mesh refinement can be restricted to specific areas that require more precise grid structure, such as the fractures network and wellbore zone (Karimi-Fard and Firoozabadi, 2003). According to the Galerkin variational method, unknown variables, in this case pore pressure, P and displacement vector, u are approximated based on the nodal values in each element.

5 × 106Pa

n

P= ∑ Pj Nj u = j =1

n

∑ ui,j Nj i = x, y, z

Where Nj is the shape function of j

Fig. 4. Numerical solution for model validation.

4

(20)

j =1 th

node. ux , u y and uz are

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

Fig. 5. Fracture distribution for Section 5.1.

Table 2 input data for Section 5.1. Rock properties Domain dimensions No-load fracture aperture, d f ,0

1 × 10−3m

Maximum fracture aperture closure, Dn

0.98 × 10−3m

{ε} = [B]{u}e

(29)

{σ} = [E ]{ϵ}

(30)

After calculating the above-mentioned matrices for each element, they are combined to form the global system of equations.

200m × 200m × 100m

Ne, f

Ne, m

[M ] =



[T ]em +d f

e =1

∑ [T ]ef e =1

Ne, f

Ne, m

Table 3 Fracture properties for different rock types.

[M ] =



[M ]em +d f

e =1

Fracture properties

Case 1

Case 2

Case 3

Initial normal stiffness, Kn, i [GPa /m] Shear stiffness, Ks [GPa /m]

1200 600

120 60

12 6

(31)

∑ [M ]ef e =1

(32)

Ne, m

[K ] =



[K ]e (33)

e =1 Ne, m

[Q ] =

displacements in x , y and z direction, respectively. Multiplying Eqs. (5) and (10) by a weight function and integrating over the domain of each element, results in Eqs. (21) and (22). In the Galerkin method this weight function is the shape function itself.

⎡⎛





∫Ω ⎢⎣α∇. ∂∂ut +S ∂∂Pt +qsc′ ⎥⎦ Nd Ω

=

[K ]{u}−[Q]{P}={ f } (21)

[M ]

[M ]

(23)

d {P} d{u} + [Q]T + [T ]{P} = {qsc} dt dt

e

(24)

Where [K ], [Q], [M ] and [T ] are coefficient matrices, that are functions of element geometry, fluid characteristics, and mechanical properties of both rock and fractures:

[K ]e =

∫Ω [B]T [E ][B] d Ω

(25)

[Q]e =

∫Ω α [B]T m [N ] d Ω

(26)

[M ]e =

∫Ω S [N ]T [N ] d Ω

(27)

[T ]e =

∫Ω kμ [∇N ]T [∇N ] d Ω

Where [m] = {1 1 1 0 0

0}T

d {P} d{u} +[Q]T +[T ]{P}={qsc} dt dt

(35) (36)

Where {P} and {u} are the global nodal values of pore pressure and displacement vector, respectively. Several approaches can be taken to solve Eqs. (35) and (36). One approach is to solve both system of equation simultaneously using a fully implicit scheme (Garipov et al., 2016). Since the equations are coupled together, this method will yield the most accurate result. However, because fracture properties are time-dependent, the coefficient matrices, must be updated in each time step. If we denote the number of total nodes by n , this means that in every time step a system of 4n must be solved simultaneously, in order to find the nodal values of pore pressure and displacement vector. Another approach is to first solve one system of equations and then solve the other set of equations based on the values acquired from the first system of equations. Since the time-dependency of coefficients has its roots in the fluid flow equation, this method will allow us to only change the coefficient matrices of the fluid flow equation. This in turn, will reduce the computational cost significantly. The method introduced in this work is similar to the latter approach. As shown by Fig. 2, in each time step, first the fluid flow equation is solved implicitly to find pore pressure, P . Then, based on the calculated values for pore pressure, P , the second set of equations is solved to find the nodal values of the displacement vector, u. After calculating the displacement vector, the stress field and then new fracture apertures are calculated using Eq. (18). Iteration is done within each time step until the unknown variables, P and u, converge.

(22)

If displacement vector and pore pressure in above equations, are replaced by the approximations of Eq. (20), they can be written in matrix form:

[K ]{u} − [Q]{P} = { f } e

(34)

Subscripts m and f denote porous media and fractures, respectively. The global system of equations is:





∫Ω ∇.⎜⎝ kμ ∇P⎟⎠ Nd Ω

[Q]e

e =1





∫Ω ⎢⎣ ⎜⎝K + 13 G⎟⎠ ∇ϵ+G∇2u−α∇P⎥⎦ Nd Ω= ∫Ω fNd Ω



(28) . Also, by definition: 5

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

Fig. 6. Steady-state pressure distribution for different types of rock.

4. Model validation

The same process is repeated for next time steps.

In this section, in order to validate the accuracy of the presented method, the numerical solution of a two-dimensional problem is compared to the solution obtained by Gu et al. (2014). First, the problem is described and then the numerical solutions are presented. The problem consists of a two-dimensional domain with two horizontal wells at each side, as shown by Fig. 3. Four fractures are located inside the domain according to Fig. 3. Injecting and producing wells, each operate at constant pressure. Pore pressure is uniform at t = 0 and in-situ stress applied to the outer boundaries of the domain is constant at all times. Also, the fractures have similar properties at t = 0 . The input data used in this section is listed in Table 1. Fig. 4 shows the numerical solution for the fracture apertures along the width of the domain in steady-state conditions. Comparison between the results obtained in the present work and the work done by Gu et al., 2014. demonstrates that the numerical method is valid and accurate. It can be seen in Fig. 4 that the variation in fracture aperture

3.2. Computational grid Acquiring accurate solutions in fractured porous media requires an appropriate grid structure. Geological features such as fractures, faults and well locations must be considered carefully, in order to create such an appropriate grid. This involves adjusting the grid resolution near such features. Optimum grid resolution can improve solution stability and accuracy (Bahrainian and Daneh Dezfuli, 2014). In the present work, Gmsh software has been used to generate three-dimensional unstructured grid (Geuzaine and Remacle, 2009). Gmsh software implements a Delaunay-based triangulation scheme for grid generation. Here, an adaptive unstructured grid, using tetrahedral elements. Since fractured porous media have complex geometries, unstructured grid is more suitable for such formations. 6

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

Fig. 7. Fracture aperture distribution for different types of rock.

introduced and the validity of its solution has been established. In this section, two problems are investigated by this method, in order to demonstrate the capabilities of the presented model.

along the width of the domain is almost a straight line with constant slope. This is due to the fact that in steady-state conditions, pressure gradient between the injecting and the producing well is approximately constant, which means that pore pressure variation is also a line with constant slope. Therefore, according to the governing equations, the change in fracture aperture would also be a line with constant slope.

5.1. Fractured reservoir with a single production well Fig. 5 shows the fracture distribution inside the reservoir. The reservoir is 200 m×200 m and has a formation thickness of 100 m. Fractures no. 1 to no. 7 are connected to each other and fracture no. 5 is a hydraulic fracture that intersects the production well. The produc-

5. Application in reservoir simulation So far, a coupled fluid flow and geomechanics model has been

Fig. 8. Left: Volumetric and Right: Cumulative production rate vs. time for different rock types.

7

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

than shear dilation effects. Fig. 8 shows the production rate and cumulative production variations with time, respectively. It is evident that approximately at day 4 fluid flow reaches steady-state and after that production rate remains constant and cumulative production is a line with constant slope. In case 1, in which fractures are most rigid, there is very little variation in production rate compared to the non-sensitive case. But in cases 2 and 3, reduction in production rate. The reason for this rate reduction is that as pore pressure of the reservoir decreases, the aperture of the fractures also decrease, making it harder for the fluid to travel through the fracture network. Since, pressure gradient between well and reservoir is similar to non-sensitive conditions, more resistance against fluid flow result in a decrease in the production rate.

Fig. 9. Schematic of reservoir in Section 5.2.

5.2. Two-layer fractured reservoir with production and injection wells

Table 4 Input data for Section 5.2.

In this case, we consider a fracture reservoir that consists of two geological layers. The lower layer is 100 times more permeable than the upper shale layer. As shown in Fig. 9, there are two horizontal wells located at opposite sides of the reservoir. The injection well is placed inside the lower layer and production well is located in the upper shale layer. Fracture network geometry and distribution is similar to case 5.1, and fluid and rock properties are given in Table 4. Outer boundaries of the reservoir are impermeable and no-flow boundary condition is imposed. Initial reservoir pore pressure is 2 × 107Pa , production well pressure is 107Pa , and injection well pressure is equal to 3 × 107Pa . Similar to Section 5.1, the outer boundaries do not move and the displacement vector is equal to zero at the outer boundaries. Fluid properties are similar to Section 4 and rock properties are given in Table 2. Other rock properties are similar to Section 4. In-situ stress is constant and equal to 5 × 107Pa . Fig. 10 displays the pressure distribution at 4 different times in a section of reservoir that intersects both wells. At t = 1hr there is still no significant amount of pressure diffusion and fluid flow has not yet reached to the fracture system. Once we get to t = 1day , it is evident from the figure that fluid is being transported via the fracture system from the injection well to the production well. At this stage, the matrix does not have a significant contribution in producing the fluid, but since the lower layer is more permeable, it plays a role in the injection process. At t = 10days the fracture system is almost completely pressurized by the injection well and its pressure is equal to the injection well. Here, it can be seen that matrix is starting to contribute in the production of fluid. When the flow reaches steady-state conditions, almost all of the reservoir has the same pressure as the injection well. The reason is that the lower more permeable layer is much more capable of transporting the fluid, but since the upper shale layer is relatively impermeable, fluid starts to build up inside the reservoir. Therefore the reservoir pore pressure also builds up until it reaches the injection well pressure. Changes in injection and production rates with time are displayed in Fig. 11. Four distinct time spans can be distinguished on the injection and production rate curves. It can be seen that since the injection well is situated in the more permeable layer, the injection rate is bigger than the production rate at all times. In the first time span (I), both rates are rapidly decreasing, due to the high pressure gradients around the wellbores at early times. In the second stage (II), production rate stays unchanged. The reason for this is that the decrease in reservoir pore pressure around the production rate is compensated by the flow coming from the injection rate and traveling through the fracture network towards the production well. Therefore, pressure gradient near the production well remains constant and as a result production rate also does not change. As flow reaches the third stage (III), the increase in pore pressure due to injection surpasses the decrease in pore pressure due to production in near production well zone. This leads to an increase in flow rate in the production well. In

Rock properties Domain dimensions Lower layer permeability, km, l

200m × 200m × 100m

5 × 10−14m2

Upper layer permeability, km, u

5 × 10−16m2

No-load fracture aperture, df , 0

1 × 10−3m

Maximum fracture aperture closure, Dn

0.98 × 10−3m

Fluid properties Viscosity, μ

1.8 × 10−3Pa. s

Initial and Boundary conditions Initial pore pressure, Pi

20 × 106Pa

Injecting well pressure, Pinj

30 × 106Pa

Producing well pressure, Pprod

10 × 106Pa

In-situ stress, σis

20 × 106Pa

tion well is positioned at (x, y)=(100m, 100m ) and is perforated from z = −25m to z = 25m . Initial reservoir pore pressure is 5 × 107Pa , the production well pressure is 107Pa , and pore pressure at the outer boundaries of the reservoir is 5 × 107Pa . It is assumed that outer boundaries of the reservoir do not move, thus the displacement vector is equal to zero at the outer boundaries. Fluid properties are similar to Section 4 and rock properties are given in Table 2. Other rock properties are similar to Section 4. In-situ stress is constant everywhere and equal to 5 × 107Pa . We consider three types of rock with three different values for normal and shear fracture stiffness, which are given in Table 3. Fig. 6 shows the comparison of steady-state pressure distribution for nonsensitive and sensitive reservoir. It can be seen that as we decrease the initial normal fracture stiffness, Kn, i and shear stiffness, Ks , pressure diffusion also decreases. This is due to the fact that with lower stiffness, fracture are more easily affected by the stress field, and therefore for lower values of Kn, i , fracture aperture reduction is more dramatic. This reduction in fracture aperture is caused by the production process. As the production process goes on, reservoir pore pressure begins to deplete, changing the stress distribution within the reservoir. This change in stress field leads to the reduction in fracture aperture. Conversely, for injection processes, fractures might open up, allowing the fluid flow to travel through the fracture system with less difficulty. Reservoir pressure depletion causes normal closure of fracture aperture and for less values of normal fracture stiffness, this closure becomes greater. However, as discussed in previous sections, shear dilation acts against this mechanism. In other words, lower shear stiffness makes it easier for the fracture planes to slide on one another. This means that for lower values of Ks , shear dilation effect is more significant. But as we see in Figs. 6 and 7, for less rigid rocks, fracture aperture closes more. The reason behind this is that, in this case, direct normal displacement of fractures due to normal stress is much greater 8

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

Fig. 10. Pressure distribution for Section 5.2.

the fourth stage (IV), fluid flow becomes steady and injection rate equals the production rate. Fig. 12 demonstrates the variations in average fracture aperture with time for every fracture. Average fracture aperture for all fractures rapidly increases until t = 100days , and after that it remains roughly the same. The increase in fracture aperture is due the pressure build up during the injection process. Since injection rate is always greater than the production rate, pore pressure of the reservoir increases, which as a result, leads to an increase in fracture aperture.

along with an adaptive unstructured grid was used to numerically solve the coupled problem. An experimental model was incorporated to relate the fracture aperture with the reservoir stress field. The accuracy of the presented model has been validated by previously published results. The method has been used to solve two problems regarding the fluid flow in fractured porous media. The following points can be inferred from the numerical results: 1. Fracture aperture is dependent on the applied stresses to the fracture surfaces. But, the degree of this dependency is highly related to mechanical properties of rock and the production/injection conditions. 2. Normal displacement of the fracture aperture due to normal stress is much more significant than normal displacement due to shear dilation effects. 3. Production rate in fracture reservoirs is very sensitive to changes in fracture aperture and permeability. Therefore, the geomechanical

6. Conclusion In this work, a coupled fluid flow and geomechanics numerical method based on the discrete-fracture model and Finite Element method is introduced. This method, utilizes a sequential implicit scheme for solving the poroelasticity and fluid flow (continuity) equations. The Galerkin method with Finite Element discretization 9

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

Fig. 11. Left: injection and production rates vs. time, right: cumulative injection and production vs. time. Gutierrez, M., Lewis, R.W., Masters, I., 2001. Petroleum reservoir simulation coupling fluid flow and geomechanics. Soc. Pet. Eng. 4, 164–172. Gu, S., Liu, Y., Chen, Z., 2014. Numerical study of dynamic fracture aperture during production of pressure-sensitive reservoirs. Int. J. Rock. Mech. Min. Sci. 70, 229–239. Hu, L., Winterfield, P.H., Fakcharoenphol, P., Wu, Y., 2013. A novel fully-coupled flow and geomechanics model in enhanced geothermal reservoirs. J. Pet. Sci. Eng. 107, 1–11. Jaffre, J., Mnejja, M., Roberts, J.E., 2011. A discrete fracture model for two-phase flow with matrix-fracture interaction. Procedia Comput. Sci. 4, 967–973. Jalali, M.R., Dusseault, M.B., 2012. Coupling geomechanics and transport in naturally fractured reservoirs. Int. J. Min. Geo-Eng. (IJMGE) 46 (2), 105–131. Karimi-Fard, M., Durlofsky, L., Aziz, K., 2004. An efficient discrete fracture model applicable for general-purpose reservoir simulators. SPE J. 9 (2), 227–236. Karimi-Fard, M., Firoozabadi, A., 2003. Numerical simulation of water injection in fractured media using the discrete-fracture model and the Galerkin method. SPE Reserv. Eval. Eng. 6 (2), 117–126. Khoei, A.R., Shamloo, A., Anahid, M., Shahim, K., 2006. The extended finite element method (X-FEM) for powder forming problems. J. Mater. Process. Technol. 177 (1– 3), 53–57. Kilmer, N.H., Morrow, N.R., Pitman, J.K., 1987. Pressure sensitivity of low permeability sandstones. J. Pet. Sci. Eng. 1 (1), 65–81. Kim, J., Moridis, G., Yang, D., Rutqvist, J., 2012. Numerical studies on two-way coupled fluid flow and geomechanics in hydrate deposits. Soc. Pet. Eng. 17 (2), 485–501. Lamb, A.R., Gorman, G.J., Elsworth, D., 2013. A fracture mapping and extended finite element scheme for coupled deformation and fluid flow in fractured porous media. Int. J. Numer. Anal. Methods Geomech. 37 (17), 2916–2936. Martin, V., Jaffre, J., Roberts, J., 2005. Modeling fractures and barriers as interfaces for flow in porous media. SIAM J. Sci. Comput. 26 (5), 1667–1691. McClure, M.W., Babazadeh, M., Shiozawa, S., Huang, J., 2016. Fully coupled hydromechanical simulation of hydraulic fracturing in 3d discrete-fracture networks. Soc. Pet. Eng. 21, 1302–1320. Minkoff, S.E., Stone, C.M., Bryant, S., Peszynska, M., Wheeler, M.F., 2003. Coupled fluid flow and geomechanical deformation modeling. J. Pet. Sci. Eng. 38, 37–56. Moinfar, A., Sepehrnoori, K., Johns, R., Varavei, A., 2013. Coupled geomechanics and flow simulation for an embedded discrete fracture model. In proceedings of: Society of Petroleum Engineers – SPE Reservoir Simulation Symposium 2, 1238–1250. Monteagudo, J.E.P., Firoozabadi, A., 2004. Control-volume method for numerical simulation of two-phase immiscible flow in two- and three-dimensional discretefractured media. Water Resour. Res. 40, 1–20. Norbeck, J., Huang, H., Podgorney, R., Horne, R., 2014. An integrated discrete fracture model for description of dynamic behavior in fractured reservoirs. In proceedings of: 39th workshop on geothermal reservoir engineering, Stanford university, California. Olsson, R., Barton, N., 2001. An improved model for hydromechanical coupling during shearing of rock joints. Int. J. Rock. Mech. Min. Sci. 38 (3), 317–329. Pedrosa, O.A., 1986. Pressure transient response in stress-sensitive formations. In: Proceedings of the SPE California Regional Meeting, Oakland, California, USA, 203214. Philip, Z.G., Jennings, J.W., Olson, J.E., Laubach, S.E., Holder, J., 2005. Modeling coupled. Soc. Pet. Eng. 8 (4), 300–309. QanbariF., ClarksonC.R., 2012. Rate transient analysis of stress-sensitive formations during transient linear flow period. In: Proceedings of the SPE Canadian Unconventional Resources Conference, Calgary, Alberta, Canada. Qanbari, F., Clarkson, C.R., 2014. Analysis of transient linear flow in stress-sensitive formations. SPE Reserv. Eval. Eng. 17 (1), 98–104. Rutqvist, J., Stephansson, O., 2003. The role of hydromechanical coupling in fractured rock engineering. Hydrogeol. J. 11 (1), 7–40. Segura, J.M., Paz, C.M., de Bayser, M., Zhang, J., Bryant, P.W., González, N.A., Rodrigues, E., Vargas, P.E., Carol, I., Lakshmikantha, M.R., Das, K.C., Sandha, S.S.,

Fig. 12. Fracture aperture vs. time.

effects on fracture network should be taken into account for accurate prediction of well performance. References Bahrainian, S.S., Daneh Dezfuli, A., 2014. A geometry-based adaptive unstructured grid generation algorithm for complex geological media. Comput. Geosci. 68, 31–37. Barton, N., Bandis, S., Bakhtar, K., 1985. Strength, deformation and conductivity coupling of rock joints. Int. J. Rock. Mech. Min. Sci. Geomech. Abstr. 22 (3), 121–140. Biot, M.A., 1941. General theory of three‐dimensional consolidation. J. Appl. Phys. 12, 155–164. Castelletto, N., White, J.A., Tchelepi, H.A., 2015. Accuracy and convergence properties of the fixed-stress iterative solution of two-way coupled poromechanics. Int. J. Numer. Anal. Methods Geomech. 39, 1593–1618. Chen, H.-Y., Teufel, L., 2000. Coupling fluid-flow and geomechanics in dual-porosity modeling of naturally fractured reservoirs – model description and comparison. In: Proceedings of SPE International Petroleum Conference and Exhibition in Mexico, Villahermosa, Mexico Chen, H.-Y., Teufel, L.W., Lee, R.L., 1995. Coupled fluid flow and geomechanics in reservoir study - i. theory and governing equations. In: proceedings of: SPE Annual Technical Conference and Exhibition, Dallas, Texas, USA Chin, L.Y., Raghavan, R., Thomas, L.K., 2000. Fully coupled analysis of well responses in stress-sensitive reservoirs. SPE Reserv. Eval. Eng. 3, 435–443. Garipov, T.T., Karimi-Fard, M., Tchelepi, H.A., 2016. Discrete fracture model for coupled flow and geomechanics. Comput. Geosci. 20 (1), 149–160. Geuzaine, C., Remacle, J.-F., 2009. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Methods Eng. 79 (11), 1309–1331.

10

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al. Cerqueira, R., Mello, U., 2016. Coupling a fluid flow simulation with a geomechanical model of a fractured reservoir. American Rock Mechanics Association. In: Proceedings of the 50th US Rock Mechanics/Geomechanics Symposium, Houston, Texas, USA, 26-29. Settari, A., Mourits, F.M., 1998. A coupled reservoir and geomechanical simulation system. Soc. Pet. Eng. 3 (3), 219–226. Shamloo, A., Azami, A.R., Khoei, A.R., 2005. Modeling of pressure-sensitive materials using a cap plasticity theory in extended finite element method. J. Mater. Process. Technol. 164–165, 1248–1257. Shaoul, J.R., Ayush, A., Park, J., Pater, C.J.D., 2015. The effect of stress sensitive permeability reduction on the evaluation of post-fracture welltests in tight gas and

unconventional reservoirs. In: proceedings of the SPE European damage conference and exhibition, Budapest, Hungary. Wang, H.-X., Wang, G., Chen, Z.-X.J., Wong, R.C.K., 2010. Deformational characteristics of rock in low permeable reservoir and their effect on permeability. J. Pet. Sci. Eng. 75 (1–2), 240–243. Zhang, M.Y., Ambastha, A.K., 1994. New insights in pressure-transient analysis for stress-sensitive reservoirs. In: Proceedings of the 69th SPE Annual Technical Conference and Exhibition, New Orleans, L.A., USA, pp. 617–628. Zienkiewicz, O., Taylor, R., 2005. The Finite Element Method for Solid and Structural Mechanics. Elsevier.

11