Abstract
The effect of geomechanics is crucial in the modeling of fractured reservoirs since fractures can be more stress-sensitive than the rock matrix. In fractured reservoirs, the microscale fractures are often homogenized into the matrix continuum leading to matrix permeability anisotropy. The embedded discrete fracture model (EDFM) is becoming increasingly popular for numerical simulation of fractured reservoirs and has been successfully employed in the simulation of coupled geomechanics and flow systems. However, the anisotropic permeability of matrix cannot be considered in traditional EDFM. An approach of coupled geomechanics and flow modeling is proposed for fractured reservoirs considering anisotropic matrix permeability. In order to calculate the fluid exchange between fractures and rock matrix in an anisotropic formation, an integrally embedded discrete fracture model (IEDFM) is used. The geomechanics model of the proposed approach uses an equivalent continuum approach, which introduces an equivalent material to represent the overall deformation of the fractured rock under normal and shear stresses. The constitutive relations of the equivalent continuum are derived rigorously from stress-strain analysis, where the stress-dependent moduli of natural fractures are included. The coupled geomechanics and flow system is solved using the fixed-stress split iterative coupling strategy with the dynamic hydraulic parameters of matrix and fractures updated separately. Several examples are performed to demonstrate the applicability of the proposed approach for modeling the coupled geomechanics and flow system in fractured reservoirs considering anisotropic permeability. The effect of anisotropy is investigated, which indicates that the dynamic behavior of a fracture is highly orientation-related due to initial stress anisotropy and matrix permeability anisotropy. Simulations also show that the anisotropic matrix permeability affects the compaction in the reservoir domain, which reflects on the performance of production.
1. Introduction
During the exploitation of hydrocarbon reservoirs, pore pressure changes caused by production can lead to rock compaction, which in return affects the hydraulic properties such as porosity and permeability. In fractured reservoirs, the effect of geomechanics is more crucial since fractures can be more stress-sensitive than the rock matrix. In order to capture the effect of geomechanics on the production of fractured reservoirs, it is necessary to couple flow and geomechanics.
The fundamentals of coupled geomechanics and flow problems date back to the concept of effective stress first proposed by Terzaghi [1]. Biot [2, 3] later extended this concept and derived the quasistatic constitutive equations that relate the strain tensor to the stress tensor and fluid pressure, which is generalized as the three-dimensional poroelastic theory. In the poroelastic theory, rock blocks are composed of compressible soil skeleton and fluid that fills the pore space. However, in a fractured porous media, the single porosity poroelastic theory is no longer applicable. According to the dual continuum concept introduced by Barenblatt et al. [4], the matrix and fractures can be envisioned as two separate continua that interact with each other. Based on the dual continuum concept, a dual porosity poroelastic theory (DPP) is proposed [5–7], where the effective stress law is modified to capture different extents of deformation of the rock matrix and fractures. Then, Wittke [8] rigorously derived the constitutive model for a fractured rock using the concept of equivalent continuum. Later, series of attempts to simulate coupled geomechanics and flow problems in fractured reservoirs using the dual continuum concept and the dual porosity poroelastic theory are performed [9–15].
Though dual continuum models provide efficient ways of modeling fractured reservoirs, these models cannot handle reservoirs with complex fracture networks composed of large-scale fractures since the real fracture geometries cannot be represented [16, 17]. The embedded discrete fracture model (EDFM) is proposed by Li and Lee [18] to improve this limitation and later successfully implemented for various fracture geometries [19–21]. In the EDFM, fractures are treated explicitly to capture the complex geometry of fracture networks. The rock matrix is discretized using Cartesian grids to keep the computational efficiency, and the fracture elements are represented by individual control volumes that are connected to the parent matrix blocks. The EDFM has become an increasingly popular method for modeling fractured reservoirs and has been successfully adopted in coupled geomechanics and flow modeling. Moinfar et al. [22] proposed a coupled flow and geomechanics model for fractured reservoirs using the EDFM. Sangnimnuan et al. [23] presented a coupled flow and geomechanics model using the EDFM to predict stress evolution in unconventional reservoirs during production. Liu et al. [24] investigated the effects of hydraulic fracture and natural fracture properties on production behavior as well as pressure and stress evolution of shale gas reservoirs using the EDFM and FEM coupled simulations. Bai et al. [25] developed a coupled compositional flow and geomechanics simulator using the EDFM and investigated the effect of confined phase behavior in shale oil reservoirs. Ren et al. [26] proposed a coupled XFEM-EDFM hybrid model for geomechanics and flow in fractured reservoirs, which uses the extended finite element method (XFEM) to model the geomechanics and the EDFM to model the flow equations. This hybrid model is later applied in the modeling of the fractured vuggy porous media [27], hydraulic fracture propagation process [28], and enhanced geothermal system [29].
Eventhough extensive research works have been done on the coupled geomechanics and flow modeling of fractured reservoirs using the EDFM, the effect of matrix permeability anisotropy is often overlooked, which might significantly affect the performance of the reservoir in coupled geomechanics and flow simulations. Naturally fractured reservoirs contain multiple scaled fractures, including a limited number of macroscale fractures and a large number of microscale fractures. Though the macroscale fractures can be treated explicitly and efficiently using the EDFM, individually treating all the large amount of microscale fractures as control volumes using the EDFM often causes large computational burden. One practical method to treat these microscale fractures is to homogenize them into the matrix continuum, which results in matrix permeability anisotropy. Traditional EDFM cannot precisely model the anisotropic flow between fracture and matrix due to the oversimplified linear pressure distribution assumption applied in the EDFM [30]. Shao et al. [31, 32] developed an integrally embedded discrete fracture model (IEDFM), where fractures are envisioned as series of point sinks to derive the semi-analytical pressure distribution around the fractures. The IEDFM allows for the calculation of matrix-fracture connection transmissibilities considering the anisotropy of rock matrix, which is used to calculate the fluid exchange between the rock matrix and fracture.
A coupled geomechanics and flow model for the simulation of fractured reservoirs considering matrix anisotropy is proposed in this study. The flow simulation uses the IEDFM to account for the fluid exchange between fracture control volumes and the rock matrix with anisotropic permeability. In the geomechanics modeling, an equivalent continuum approach is used to derive the stiffness matrix required in the finite element method, where an equivalent material is introduced to represent the overall deformation of fractured rock. An empirical joint model is employed to capture the stress-dependent fracture moduli. Several examples are performed to demonstrate the applicability of the proposed approach for coupled geomechanics and flow modeling of fractured reservoirs with anisotropic permeability, and the effect of anisotropy in coupled geomechanics and flow systems is investigated.
2. Governing Equations
The governing equations of a coupled geomechanics and flow system consist of two parts: the mass conservation equation and the force balance equation.
2.1. Mass Conservation Equations
The mass conservation equations in fractured media are given bywhere represents the fluid phases, is the porosity, is the saturation of phase , is the density of phase , is the velocity of phase , and is the source term of phase . This mass balance equation applies to both matrix and fracture control volumes.
The velocity of phase is defined by Darcy’s law:where is the absolute permeability tensor, is the relative permeability of phase , is the viscosity of phase , is the pressure of phase , is the gravitational acceleration, and is the depth.
Discretizing equation (1) using the finite volume method in space and first-order scheme in time giveswhere subscript denotes the values of element , superscript represents the current time level, superscript represents the previous time level, is the flow term for phase between element and element , is a set that contains all the neighbor elements connected to element , and is the source term for phase .
The fluid flow through matrix-matrix, fracture-fracture, and matrix-fracture connections is universally described by the flow term , which can be expressed aswhere is the transmissibility of the connection, which is the harmonic average of the half-transmissibility of element and element :where and are the half-transmissibilities, and are the absolute permeabilities along each element center to the interface, and are the distances from the centers to the interface, and is the interface area.
2.2. Force Balance Equations
The vectors and tensors are expressed in components. The force balance equation can be obtained from linear momentum conservation law, which giveswhere is the density of the fluid-solid mixture, is the density of rock skeleton, is the acceleration vector of the rock skeleton, is the total stress tensor, and is the body force vector.
In a fractured porous media, the rock matrix and fracture are envisioned as two separate continua for the flow system; thus, the traditional effective stress law is not appropriate. The effective stress law for fractured media can be derived by extending the single porosity effective stress law, which giveswhere is the effective stress tensor, is the average pressure of matrix blocks, is the average pressure of fractures, is the Kronecker delta, is the Biot coefficient of matrix, is the Biot coefficient of fracture, is the bulk modulus of porous media, is the modulus of solid grain, and is the drained bulk modulus of porous media without fractures.
Since the loadings caused by fluid depletion do not change drastically, deformation in the reservoir can be assumed to be quasistatic. In the case of static analysis, the dynamic term in equation (7) can be subtracted. Then, substituting equation (8) into equation (7) gives
In the case of an elastic tensor, the constitutive relations are given bywhere is the fourth-order modulus tensor, and is the strain tensor.
Spatially discretizing equation (9) in incremental form using the finite element method gives the discretized force balance equation:where is the nodal stiffness matrix, is the strain-nodal displacement matrix, is the elastic stiffness matrix, which is the matrix form of the fourth-order modulus tensor , is the nodal displacement increment vector, is the volumetric stiffness vector, and are the pressure increments in matrix and fracture, and is the external loading vector.
3. Modeling of Flow in Fractured Reservoirs considering Matrix Permeability Anisotropy
The IEDFM is adopted to model the fluid flow in fractured reservoirs considering matrix permeability anisotropy. In the IEDFM, rock matrix is discretized by Cartesian grids, and additional elements are added to represent the fracture control volumes. There are three kinds of connections between control volumes in the IEDFM: matrix-matrix (M-M) connections, fracture-fracture (F-F) connections, and matrix-fracture (M-F) connections.
For M-M connections, the multipoint flux approximation (MPFA) method [33] is used to establish the connections between matrix elements in an anisotropic reservoir. Since the permeability difference between two connected matrix elements is relatively small, the pressure distribution inside matrix elements can be assumed to be linear. Thus, the parameters in equation (6) have clear physical meanings: and are the distances from the matrix element centers to the interface, and and are the matrix permeabilities that can be calculated from the absolute permeability tensor using MPFA.
For F-F connections, the linear pressure distribution is also appropriate since the permeability difference is small. The parameters are as follows: and are the distances from the fracture element centers to the common interface, and and are the fracture permeabilities .
For M-F connections, due to the anisotropy of matrix permeability, simplified linear assumptions cannot be used for the pressure distribution around a fracture. In the IEDFM, the M-F transmissibilities are calculated semi-analytically to address the anisotropic matrix permeability. A fracture is considered as a series of point sinks , as shown in Figure 1. In a 2D anisotropic medium, if a point sink locates at the origin , the pressure at position is given bywhere is the point sink flow rate, is the strata height, is the equivalent permeability of the anisotropic matrix, is the equivalent distance, and is a constant number.
In a matrix element, the pressure of point can be derived as the superposition of all the pressure drops caused by point sinks:where is the number of point sinks, is the flow rate of sink , is the equivalent distance between point and sink , and is a constant related to the fracture pressure.
The fractures are assumed to be equipotential; thus, several points on the fracture surface can be selected as the reference points, and these points have the same pressure . Applying equation (15) on every reference point forms a linear equation system:where is the equivalent distance between reference point and sink .
The linear equation system can be rewritten aswhere is defined by
By solving the linear equation system, the exact expression of pressure at point can be obtained:where is the matrix pressure which is defined as the average of all point pressures, and is the matrix volume.
The transmissibility between the fracture element and the neighboring matrix element can be derived using the matrix pressure, fracture pressure, and the flow term definition given in equation (4), which gives
The transmissibility given in equation (22) is calculated semi-analytically, which improves the accuracy and allows for anisotropic matrix permeability. Validation and detailed discussion of the IEDFM in anisotropic formations can be found in the work of Shao et al. [32].
4. Modeling of Fractured Rock Deformation
4.1. Equivalent Continuum Model
To derive the nodal stiffness matrix in equation (11), the elastic stiffness matrix is required, which describes the deformation of a fractured rock.
A fractured rock consists of two separate continua, rock matrix and fracture, which behave differently under stress. The concept of equivalent continuum is introduced to establish the stiffness matrix for the fractured rock. In the equivalent continuum model, an equivalent material is introduced to represent the overall deformation of a fractured rock. The deformation behavior of the introduced equivalent material is the combination of rock matrix deformation and the displacement of fracture planes.
The constitutive relations of the equivalent material can be rigorously derived from stress-strain analysis. In case of a fractured rock crossed by a single fracture, a coordinate transformation is first performed to simplify the derivation. After the coordinate transformation, the Z-axis of the coordinate system is set to align with the normal direction of the surface of the fracture. All kinds of normal and shear loadings are applied on the rock block to obtain the stress and strain tensors, as shown in Figure 2. Young’s modulus and shear modulus of the fracture are assumed to be much smaller than the rock matrix, and the volume fraction of fracture in a fractured rock is much smaller than the rock matrix. The rock matrix is assumed to be elastic isotropic. In the local coordinate system of fracture, the compliance matrix can be written aswhere is the volume fraction of the rock matrix in the fractured rock, is the volume fraction of the fracture in the fractured rock, is Young’s modulus of rock matrix, is the shear modulus of rock matrix, is Poisson’s ratio of rock matrix, is the compression modulus of the fracture, and is the shear modulus of the fracture. The detailed derivation of equation (23) from stress-strain analysis is given in the Appendix.
The compression modulus describes the modulus of the laterally confined fracture medium, which is given bywhere is Young’s modulus of the fracture and is Poisson’s ratio of the fracture.
Since the volume fraction of the rock matrix is much larger than the fracture, it can be assumed that . The compliance matrix of the equivalent material can be rewritten aswhere and are the compliance matrixes for the rock matrix and the fracture given by
The compliance matrix given in equation (25) is derived in the local coordinate system where Z-axis aligns with the normal direction of the fracture. To obtain the compliance matrix in the global coordinate system , a coordinate transformation needs to be performed:where is the coordinate transformation matrix from the local coordinate system to the global coordinate system, which is determined by the orientation and dip angle of the fracture.
Substituting equation (25) into equation (28) giveswhere remains unchanged under the coordinate transformation due to the assumption of elastic isotropic rock matrix.
When a fractured rock is crossed by more than one fractures, the compliance matrix of the equivalent material is derived by sequentially accumulating the contribution of each individual fracture. The overall compliance matrix of a fractured rock containing multiple fractures is given bywhere is the number of fractures, is the volume fraction of fracture , is the compliance matrix of fracture , and is the coordinate transformation matrix from the local coordinate system of fracture to the global coordinate system.
The elastic stiffness matrix in equation (11) is the inversion of the overall compliance matrix :
4.2. Constitutive Model of Fractures
The deformation of fractures under stress is significantly different from intact rock matrix, which is more complex and strongly nonlinear. Thus, a fracture constitutive model is needed to calculate the dynamic fracture modulus used in equation (27). For natural fractures, an empirical model can be employed, which is based on laboratory experimental data on interlocked rock samples with natural unfilled joints [34].
For a fracture under normal stress, Young’s modulus is given bywhere is the initial fracture normal stiffness, is the fracture aperture, is the effective normal stress, and is the maximum normal closure of fracture.
The initial fracture normal stiffness can be obtained from the following equation:where JRC is the joint roughness coefficient, JCS is the joint compressive strength, and is the zero-stress fracture aperture. In this work, the term “zero-stress” refers to the state when no stress is applied, and the term “initial” refers to the state under the initial stress field before the simulation begins.
The maximum normal closure of fracture can be estimated from the following empirical relation:
The fracture aperture is stress-dependent, and the fracture aperture change under normal stress is given by
Fracture behavior under shear stress is more sophisticated and less well understood compared to deformation under normal stress. For simplicity, a nonlinear model that only applies to the prepeak fracture shear deformation is used to describe fracture behavior. The shear modulus of fracture is obtained from the empirical relation:where is the stiffness number, is the stiffness exponent, is the shear stress, is the peak shear stress, and is the failure ratio.
The stiffness number can be estimated using the empirical relation:
The peak shear stress is given bywhere is the fracture residue angle.
4.3. Validation of the Fractured Rock Geomechanics Model
The proposed geomechanics model of the fractured rock is compared with the experimental data from the work of Bandis et al. [34]. The experimental data are acquired from applying gradually increasing loading on the top of single-jointed rectangular rock samples. The curves are drawn from experiments on a limestone sample with a fully interlocked joint in the first loading cycle. Young’s modulus for limestone matrix is 49.0 GPa. The joint roughness coefficient (JRC) is 7.6, and the joint compressive strength (JCS) is 157 MPa. The zero-stress aperture of the fracture is 0.2 mm.
Figure 3 shows the comparison of the normal deformation of solid rock, fracture, and fractured rock between experimental data and calculated results. The good agreement between experimental and calculated data validates the equivalent continuum model and the fracture constitutive model. The deformation of solid rock is approximately linear, indicating the rock matrix can be assumed to be elastic. The deformation of the fracture shows nonlinear behavior and almost stops increasing at high stress level, since the maximum normal closure is about to be reached. The deformation of fractured rock is a combination of fracture and solid rock deformation, which verifies the assumption of the equivalent continuum model.
5. Coupling Strategy and Parameter Update
The coupled geomechanics and flow problem is discretized using a hybrid finite element-finite volume method (FE-FVM). The FVM is applied for the flow model with the primary variables located at the grid center, which is given in equations (3) and (4). The FEM is applied for the geomechanical model with the displacement vectors located at the grid nodes, which is given in equation (11). The coupled geomechanics and flow problem is solved using the iterative coupling strategy.
5.1. Fixed-Stress Split Iterative Coupling Strategy
In the iterative coupling strategy, the coupled problem is split into flow and mechanics subproblems; then, flow and geomechanics modules can be invoked sequentially to solve the coupled problem. At each iteration step, either flow or mechanics problem is solved first, and then, the other problem is solved using the solution information of the previous problem. This sequential process is repeated until the solution meets the convergence criteria, and the simulation moves on to the next time step.
There are several schemes to develop an iterative coupling simulator based on the operator splitting strategies for the coupled flow and mechanics problem. The fixed-stress split scheme is employed in this work, which is recommended by Kim et al. [35]. In the fixed-stress split scheme, the flow problem is solved first while freezing the total mean stress. This scheme is unconditionally stable, and the convergence is proved by many previous works [36, 37]. This scheme is successfully applied in various attempts of solving coupled flow and mechanics problems [38, 39]. Using the fixed-stress split iterative coupling strategy, a flowchart showing the framework of solving the coupled problem is shown in Figure 4.
5.2. Parameter Update after an Iteration Step
In every iteration step, the flow simulation is first performed, and then, the geomechanics module is called after the flow simulation converges to calculate changes in stress and strain fields using the information from the flow problem solution. According to the fixed-stress split scheme, changes in matrix parameters, such as porosity and permeability, are updated by freezing the total mean stress before the next iterative step.
The variation of true matrix porosity in a deformable porous media can be derived as [40]where is the mean total stress.
By freezing the mean total stress, the reservoir porosity of rock matrix can be updated by [41]where the superscript denotes the previous iteration step, the superscript denotes the current iteration step, is the volume strain, and is the reservoir porosity, which is the ratio of the pore volume in the deformed configuration to the total volume in the undeformed configuration.
The matrix permeability tensor is updated using an empirical relation [42]:where the superscript 0 denotes the zero-stress condition and is a parameter for permeability update.
5.3. Parameter Update after a Time Step
After the solution meets the convergence criteria, the hydraulic parameters such as porosity and permeability for matrix and fractures are updated for the next time step. Since the fracture moduli are stress-dependent, the fracture stiffness is also updated to build the equivalent constitutive matrix in the next time step.
According to equation (39), the updation of matrix reservoir porosity is given bywhere the superscript denotes the previous iteration step, and the superscript denotes the current iteration step.
The matrix permeability tensor is updated by
Normal effective stress on a fracture is calculated from the stress field and fracture orientation and then used to calculate fracture aperture change using equation (35). Porosity and permeability of fractures are updated using the updated fracture apertures:
Using the updated matrix permeability tensor, the M-F transmissibilities are updated by
The fracture moduli are updated using equations (32) and (36) before the next time loop.
5.4. Validation of Coupled Geomechanics and Flow Model
The accuracy of the discretized coupled geomechanics and flow model is examined by the 1D Terzaghi’s consolidation problem [1].
Terzaghi’s problem consists of a soil column fully saturated with water, as shown in Figure 5. A uniformly distributed load is applied on the top surface of the soil column. Drainage is only allowed on the upper boundary of the column. The lateral boundaries are fixed on horizontal direction, and the bottom boundary is fixed on both horizontal and vertical directions. In this case, the computational domain is 5 m in length and 10 m in height. The parameters used in Terzaghi’s problem are given in Table 1. The grids are shown in Figure 6.
The calculated results of pore pressure and vertical displacement distributions along the central column of grids at several times are compared with the analytical solutions [43], as shown in Figure 7. A great consistency between the analytical solutions and numerical solutions can be found, which validates the proposed coupled geomechanics and flow model.
(a)
(b)
6. Example Simulations
A coupled geomechanics and flow simulator is developed to demonstrate the feasibility of the proposed model. The main programs of the flow module and mechanics module are developed on FORTRAN platform. Several examples are performed to investigate the effect of reservoir anisotropy in the modeling of the coupled geomechanics and flow system in fractured reservoirs.
6.1. Example 1: Effect of Anisotropy on the Closure of Fractures
During the depletion, the dynamic behaviors of fractures are largely orientation-related. The initial stress field in the reservoir can be anisotropic since it is affected by the geological conditions of the reservoir. Besides, the matrix permeability of the reservoir can be anisotropic due to the microscale fractures that are homogenized into the matrix continuum. To investigate these effects, a simple fractured reservoir is designed with three natural fractures, as is shown in Figure 8.
The reservoir is 420 m in length, 420 m in width, and 40 m in depth. The reservoir domain is discretized by Cartesian grids, as shown in Figure 9. The three fractures have the same length and are equally distant from the production well. Fracture 1 is parallel to the X-axis, fracture 2 is parallel to Y-axis, and fracture 3 is 45° inclined. The parameters of fractures are given in Table 1. The zero-stress aperture of fractures is 0.2 mm.
The porosity in the reservoir is 0.30. The reservoir is saturated with oil with the density of 800.0 kg/m3 and viscosity of 10.0 cP. The geomechanical configuration is shown in Figure 10. The displacement boundary conditions are as follows: the lateral boundaries are fixed on horizontal direction and the bottom boundary is fixed on all directions. The mechanical parameters are given in Table 2. The initial pressure in the reservoir is 10.0 MPa. A production well located at the center of the reservoir drains at fixed bottom hole pressure of 1.0 MPa. Constant pressure boundary condition is considered in this case; thus, the oil pressure at the boundary of the reservoir remains 10.0 MPa throughout the simulation process. The total simulation time is 3000 days.
The simulation is first run using isotropic initial stress field and matrix permeability as control. The initial stress is assumed to be 4.0 MPa, and the zero-stress permeability is 0.1 D. Under the initial stress, the initial matrix permeability is 0.088 D. Figure 11 shows the fracture aperture evolution histories of the three fractures. The curves are in perfect agreement, since the distances from the fracture centers to the production well are the same.
6.1.1. Effect of Initial Stress Anisotropy
The effect of initial stress anisotropy is investigated. Figure 12 shows the two anisotropic initial stress fields considered in this part, and the matrix permeability is kept isotropic.
Figure 13 shows the fracture aperture histories of three fractures in the two anisotropic cases. Due to the existence of the initial stress, the initial aperture of the three fractures are smaller than the zero-stress fracture aperture. For example, when , the largest closure occurs on fracture 1, which is perpendicular to the maximum initial stress direction. Fracture 1 remains only 52.6% of the zero-stress aperture before the simulation. Fracture 2, which is perpendicular to minimum initial stress direction, undergoes the smallest closure and remains 64.3% of the zero-stress aperture. Fracture 3 remains 56.2% of the zero-stress aperture. At the end of the simulation, fracture 1 remains 51.7% of the zero-stress aperture, fracture 2 remains 60.0% of the zero-stress aperture, and fracture 3 remains 54.4% of the zero-stress aperture. Because of the symmetry in geometry, the fractures behave oppositely when .
(a)
(b)
Therefore, fractures perpendicular to minimum initial stress direction have larger permeabilities before the simulation, yet are more prone to closure during the depletion, such as fracture 2 when and fracture 1 when .
6.1.2. Effect of Matrix Permeability Anisotropy
Similarly, the effect of matrix permeability anisotropy is investigated. The initial stress field is assumed to be isotropic with the value of 4.0 MPa. The two anisotropic permeability tensors under initial stress are shown in Figure 14.
Figure 15 shows the fracture aperture histories of three fractures using two anisotropic permeability tensors. Since the isotropic initial stress tensor is used, the initial apertures of the three fractures are the same. It can be found that the closure of fractures perpendicular to the maximum permeability direction is the largest, and the closure of fractures perpendicular to the minimum permeability direction is the smallest.
(a)
(b)
6.1.3. The Combined Effect of Initial Stress Anisotropy and Matrix Permeability Anisotropy
Figure 16 shows the comparison of the aperture evolution of fracture 1 with all the combinations of initial stress anisotropy and matrix permeability anisotropy discussed in this section. For ease of comparison, the ratio of fracture aperture to the initial aperture, , is used to represent the closure of fracture 1.
As shown in Figure 16, it can be concluded that a fracture undergoes the largest closure when it is perpendicular to the minimum initial stress direction, and the minimum initial stress direction aligns with the maximum permeability direction. Conversely, when a fracture is perpendicular to the maximum initial stress direction and the maximum initial stress direction aligns with the minimum permeability direction, the fracture closure is the smallest.
6.2. Example 2: A Synthetic Naturally Fractured Reservoir
Another 45 natural fractures with random length and orientations are added in the reservoir of Example 1, as shown in Figure 17. The fractures are assumed to be vertical and may only penetrate specific layers of the reservoir. The lower-middle and upper-middle layers are assumed to be penetrated by more fractures than the bottom and top layers.
With the additional fractures, the production well is now perforated on the fracture elements in the center of each layer, and the equivalent permeability of the reservoir is reduced to 0.1 mD. For the comparison of matrix permeability anisotropy, the initial stress field is fixed with x-direction initial stress at 6.0 MPa and y-direction initial stress at 2.0 MPa.
Figures 18–20 show the aperture change histories of fracture 1, fracture 2, and fracture 3, which indicate that the conclusions in Example 1 still hold for more complex fracture geometry. The closure of fractures perpendicular to the maximum permeability direction is the largest, and the closure of fractures perpendicular to the minimum permeability direction is the smallest. The aperture change of fracture 3 is only modestly influenced by the anisotropic permeability tensor since it aligns with the angular bisector of the principle permeability tensor directions.
The oil pressure profiles using three permeability tensors are shown in Figure 21. It can be found that permeability anisotropy affects the pressure distribution in the reservoir. In the direction of the maximum permeability tensor principle axis, the pressure drop is the fastest, which induces the largest aperture change in the fractures that are perpendicular to this direction.
(a)
(b)
(c)
Since the oil pressures are used to calculate the stress, strain, and displacement in the reservoir, the permeability anisotropy directly influences the results of geomechanics simulation. The distributions of mean effective stress in the reservoir using the three permeability tensors are compared, as shown in Figure 22. Stress is larger in the vicinity of fractures, which can be attributed to both the larger pressure drop around the fractures and the stiffness reduction of a fractured rock compared to intact rock.
(a)
(b)
(c)
The influence of matrix permeability anisotropy on fracture deformation and the geomechanics behavior of the reservoir is eventually reflected on the performance of the production well. Figure 23 shows the cumulative oil produced using three permeability tensors. For comparison, the results of uncoupled simulations are also included in the figure. At 1000 days, the cumulative production is decreased by 24.4%, 25.0%, and 25.3% for , , and in coupled simulations.
6.3. Example 3: Water Flooding in a Five-Spot Well System
A five-spot well system water flooding case is considered. Figure 24 shows the configuration of a synthetic 2D reservoir. The primary natural fractures (marked in red) are modified from reported seismic survey data [44], and the secondary natural fractures (marked in blue) have random lengths and orientations. The parameters of primary and secondary fractures are given in Table 3. The zero-stress aperture of primary and secondary fractures are 0.3 mm and 0.15 mm. All the fractures are assumed to penetrate the entire depth of the reservoir.
(a)
(b)
The reservoir is 500 m in length, 500 m in width, and 40 m in depth. The porosity of the reservoir is 0.20, and the initial oil saturation in the reservoir is 0.8. The initial pressure in the reservoir is 15.0 MPa. In this example, an isotropic initial stress field is applied. The initial stress is 4.0 MPa in the x-direction and y-direction and 19.0 MPa in the z-direction.
Four cases with different permeability tensors are investigated. The zero-stress permeabilities tensors for four cases are given in Table 4. The four tensors have the same equivalent permeability of 10.0 mD. Under the initial stress, the components of permeability tensors are reduced, and the equivalent permeability is reduced to 8.23 mD. Figure 25 shows the permeability tensors under initial stress.
Water is injected from an injection well in the center of the reservoir at the fixed pressure of 15.0 MPa. Four wells located at four corners of the reservoir produce the fixed bottom hole pressure of 5.0 MPa. The water density is 1000.0 kg/m3, and the oil density is 800.0 kg/m3. The water viscosity is 0.9 cP, and the oil viscosity is 3.0 cP. The relative permeability curves are assumed to be linear. A closed boundary is considered in this example. Therefore, all the boundaries of the reservoir are impermeable. The geomechanical configuration in this example is the same with that of Example 1. The simulation is run for 10000 days.
Figure 26 shows the oil saturation profiles for four anisotropic cases after 10000 days of simulation, which clearly demonstrate the effect of matrix permeability anisotropic on the water displacement process. A two-dimensional permeability tensor has two principle axes, and the direction of the maximum permeability tensor principle axis becomes a favored direction for water flooding in anisotropic reservoirs.
(a)
(b)
(c)
(d)
The saturation distribution shows that the existence of fractures causes strong heterogeneity across the reservoir. Due to the aperture difference between primary and secondary fractures, the primary fractures have larger impacts on the oil saturation distribution. The primary fractures become major flow channels in the reservoir, and the secondary fractures increase the heterogeneity of the reservoir. Fractures whose orientation aligns with the direction of the maximum permeability tensor principle axis create fast water displacement channels.
Figure 27 shows the comparison of the oil and water producing rate of four cases in four production wells. The differences in producing rates result both from the anisotropic matrix permeability tensors and the heterogeneity induced by fractures. Production wells in the direction of the maximum permeability tensor principle axis to the injection well have the highest producing rate, such as well 1 and well 3 in Case 3 and well 2 and well 4 in Case 4. Production wells in the direction of the maximum permeability tensor principle axis to the injection well have the highest producing rate, and those in the direction of the minimum permeability tensor principle axis to the injection well have the lowest producing rate.
(a)
(b)
(c)
(d)
Figure 28 shows the oil pressure distribution after 10000 days. Pressure drop in the direction of maximum permeability tensor principle axis is the slowest, and pressure drop in the direction of minimum permeability tensor principle axis is the fastest. The symmetry of the pressure fields is broken by the fractures. The geometry of fractures also influences the overall pressure depletion in the reservoir since the preferential flow path created by fractures induces larger water inflow to compensate for the pressure drop. The average oil pressures in the reservoir are 10.19 MPa for Case 1, 10.13 MPa for Case 2, 10.28 MPa for Case 3, and 10.38 MPa for Case 4.
(a)
(b)
(c)
(d)
Due to the differences in pressure profiles, the deformation distributions in the reservoir in four cases are also different, which is illustrated in the volume strain profiles shown in Figure 29. The average volume strains are 6.41E – 4 for Case 1, 6.49E – 4 for Case 2, 6.30E – 4 for Case 3, and 6.17E − 4 for Case 4.
(a)
(b)
(c)
(d)
7. Conclusions
A coupled geomechanics and flow simulation approach is proposed for modeling fractured reservoirs considering matrix permeability anisotropy. The connections between matrix elements are established using the multipoint flux approximation (MPFA) to adapt the matrix permeability anisotropy. Simulation of the fluid exchange between the matrix and fractures uses the integrally embedded discrete fracture model (IEDFM) to account for the anisotropic permeability of the rock matrix. The geomechanical model of the approach uses an equivalent continuum approach. The compliance matrix is derived through stress-strain analysis, and the stress-dependent moduli of fractures are derived using an empirical model. In the proposed approach, flow simulation and geomechanics simulation are performed sequentially using the fixed-stress split iterative coupling strategy. In the iteration loops, the matrix porosity and permeability are updated using the results of pore pressure and strain field. The fracture porosity and permeability are updated using the aperture change calculated from normal stress.
Several examples are presented to demonstrate the applicability of the proposed approach for simulating the coupled geomechanics and flow system in fractured reservoirs with anisotropic rock matrix. The effect of reservoir anisotropy in coupled geomechanics and flow modeling is investigated. The numerical results indicate that the dynamic behaviors of fractures are largely orientation-related due to the initial stress anisotropy and matrix permeability anisotropy. Fractures perpendicular to the minimum initial stress direction are more fluid conductive yet more prone to deformation due to depletion, and fractures perpendicular to the maximum matrix permeability direction undergo the largest closure due to the rapid pressure drop. The distributions of oil pressure and saturation show that the presence of fractures causes strong heterogeneity in fractured reservoirs. Fractures that align with the direction of the maximum permeability tensor principle axis create fast water displacement channels in the reservoir. Besides, the matrix permeability anisotropy affects the compaction in the reservoir domain, which eventually reflects on the performance of production.
Appendix
As shown in Figure 30, the stress-strain analysis process for the equivalent material is illustrated. With normal effective stress applied on the top and bottom surface, the fractured rock undergoes both vertical and lateral deformations. Due to the geometric symmetry of the equivalent material, the x-direction deformation is considered as a representative of the lateral deformation.
The normal strain in the z-direction of the equivalent material is the volume average of the normal strain of rock matrix and fracture:where is the normal strain of the equivalent material, is the normal strain of the rock matrix, and is the normal strain of the fracture.
Since the volume fraction of the rock matrix is much larger than the fracture, the lateral strain can be assumed to be equal to the lateral strain of the rock matrix, and the fracture is considered to be laterally confined in the rock matrix. The relation between z-direction normal effective stress and x-direction normal strain can be given using the elastic constitutive model of rock matrix:where is the normal effective stress of the equivalent material.
Using the assumption of complete fluid exchange between rock matrix and fracture, the normal effective stress in the z-direction is equal for rock matrix and fracture:where is the normal effective stress of the rock matrix and is the normal effective stress of the fracture.
Equations (A.1) and (A.3) combined with the elastic constitutive models of the rock matrix and fracture further give the relation between the z-direction normal effective stress and the z-direction normal strain:where the modulus of fracture uses the compression modulus, which describes the normal deformation of a lateral confined material.
Equations (A.2) and (A.4) fully describe the stress-strain relations of the equivalent material under z-direction normal effective stress. Similar derivation can be applied for all kinds of normal and shear stresses. The components of the equivalent material strain tensor can be given bywhere is the shear strain of the equivalent material, is the shear strain of the rock matrix, and is the shear strain of the fracture.
The components of the equivalent material stress tensor can be given bywhere is the shear stress of the equivalent material, is the shear stress of the rock matrix, and is the shear stress of the fracture.
The constitutive relations of the equivalent material can be derived from the stress and strain tensors, which is rewritten as the compliance matrix in equation (23).
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Acknowledgments
This project was funded by the Joint Fund for Enterprise Innovation and Development of NSFC (U19B6003-02-06) and National Natural Science Foundation of China (51674010).