# Study on Stress-Dependent Permeability of Fracture Networks in Fractured Porous Media

Authors: Yueli Feng; Yuetian Liu; Gang Lei
Journal: Geofluids
 (2021)

---

## Abstract

In order to investigate the stress-sensitive characteristics of fracture networks under reservoir actual stress condition and its influence on the seepage in fractured porous media, we carried out permeability tests on experimental models with fracture networks under constant-volume boundary condition. In addition, a novel analytical stress-dependent permeability model of fracture networks in different directions was derived. Based on the test results and the proposed analytical model, the effects of various parameters (e.g., initial fracture aperture, fluid pressure, rock elastic modulus, effective-stress coefficient, and fracture dip) on deformation characteristics of fracture networks and the corresponding permeability tensor of fracture networks were studied. The research results show that, for a fractured porous media with a single group of fractures, the principal value of permeability is always parallel to the fracture-development direction. With increasing effective stress, the principal value of permeability decreases; however, the principal value direction remains unchanged. Moreover, for the fractured porous media with multiple sets of fractures, the principal direction of equivalent permeability will be inclined to the fractures with larger fracture aperture. Specifically, for the fractured porous media with two sets of intersecting fractures, the principal direction of equivalent permeability is parallel to the angular bisector of these two sets of intersecting fractures. Furthermore, the greater the difference of the fracture aperture change rate under effective stress, the more obvious the deviation of the permeability principal direction. The derived analytical model is of great theoretical and scientific significance to deepen the understanding of the stress-sensitive permeability of fractured reservoirs.


---

## Body

## 1. Introduction

When reservoir pressure changes, fractures in different directions will show different deformation characteristics, resulting in a complex oil-water seepage law, which seriously affects the development effect of water drive in fractured reservoirs. The research on the influence of fracture direction on fracture pressure-sensitive deformation characteristics and seepage law of fractured reservoirs can provide a new idea for the development adjustment of similar reservoirs.It is well acknowledged that fractured rock media generally have strong stress sensitivity, and the stress sensitivity characteristics of fractured porous media vary with different constraints. Physically, under actual reservoir condition, the fractured rock medium may be subjected to multiple stress and strain boundary conditions [1–5]. In general, the fixed-stress boundary condition and the constant-volume (fixed displacement) boundary condition are the most important two boundary conditions, which have been widely studied in scientific and engineering fields, such as physics, hydraulics, chemistry, petroleum, and engineering [6–13]. Physically, during the process of reservoir development, with an increase in effective stress, physical properties (e.g., porosity and permeability) of the rock will change. Cui et al. [14] concluded that, due to increasing of effective stress, rock permeability decreases sharply. In addition, they commented that, under a given effective stress, permeability loss of fractured porous media was more serious than that of conventional reservoirs, which could not be ignored. Since permeability is one of the key factors for characterizing reservoir production, understanding the stress-dependent permeability of fractured porous media is helpful for reservoir dynamic analysis and development plan optimization.Currently, many scholars [6–9, 14, 15] have carried out a large number of stress sensitivity experimental studies on fractured rock media under constant-stress boundary conditions. Based on these physical experimental studies, many mechanisms about the interaction between porous media and fluids have been investigated, and some theoretical understandings have been achieved. For example, Palmer and Mansoori [15] obtained a classical model to calculate permeability, which provides a function for effective stress and matrix shrinkage, and the equation is suitable for uniaxial-strain condition. Wang et al. [6] used physical simulation results to find the relationship between permeability evolution, applied stress boundaries, pore pressure, and deviatoric stress. Palmer [7], Shi and Durucan [8], and Cui and Bustin [9] determined the compressibility deformation for a fracture system and a matrix block, which follows the change of effective stress under constant-stress boundary. Cui et al. [14] considered the effect of adsorption-induced volumetric strain in the models related to confining pressure and permeability. However, successful reports of stress sensitivity experiments of fractured media with constant-volume boundary conditions are rare.Constant-stress boundary conditions that are equivalent to the constant total stress apply to reservoir boundaries. When fluid pressure in the pores changes, the effective stress is equal to the difference between total stress and pore-fluid pressure [16]. In fact, during the process of reservoir development, when pore-fluid pressure decreases, both the reservoir and surrounding rock (overlying strata, nonreservoir outside the outer boundary of a reservoir, and underlying strata) will not be able to maintain the initial stress balance. As a result, the internal stress of the reservoir will redistribute until it reaches a new stress balance state [17]. Therefore, the constant-stress condition cannot fully reflect the deformation characteristics of a pressure-sensitive reservoir. When the reservoir boundary condition has a constant volume (zero displacement boundary) and the pore-fluid pressure changes, both the effective stress and the total stress for the reservoir boundary change. The pressure-sensitive deformation characteristics of fractures, which were obtained under this boundary condition are much closer to the actual reservoir. Therefore, the characteristics of reservoir pressure sensitivity for a constant-volume boundary condition are required to clearly understand the directional pressure-sensitivity effect of fractures. This is needed to further improve the permeability characterization of fractured media. Lu et al. [1] stated that the model (for the constant-volume condition) represents a special case of the constant effective-stress model.To confirm the theory of pressure sensitivity in fractured media, Massarotto et al. [18], Ma et al. [2], Wang et al. [3], Ma et al. [19], and Lu et al. [1] used the constant-volume condition to study fractured media deformation. Massarotto et al. [18] discussed the geomechanical theories for reservoir stiffness and the pressure-arching effect. The group used coal-mine proofs to provide support for the development of important principles related to the “initial volume” condition of deep fractured coal seams, rather than constant external pressure conditions. Ma et al. [2] proposed a new model, which was based on the principle of volume balance. The group used the volume theory to calculate the volume balance between fractured coal media, solid particles, and pores. The advantage of the proposed model was that all calculation input-parameters were easy to measure. It was assumed that the constant-volume condition satisfies the uniaxial-strain condition under the condition that the vertical strain was zero. In addition, Ma et al. [2] believed that the question whether the fractured reservoir was under uniaxial strain or had a constant-volume has not truly been resolved. Inconsistent results were published, such as the research of Shi and Durucan [20] which contradicted the research of Massarotto et al. [18]. Therefore, Ma et al. [2] believe that efforts should be made to verify whether the deformation analysis of fractured media actually meets the standard of constant-volume conditions.In the studies of most researchers, fractures were rarely studied as independent objects [3, 19, 21–25]. Therefore, it is difficult for us to fully understand the characteristics of stress-sensitive permeability changes in fractured reservoirs [26].This paper is aimed at clarifying both the deformation characteristics and mathematical representation model for an independent fracture system and constant-volume boundary conditions. We used the physical simulation test results to verify the reliability of the mathematical model. Finally, the effects of different fracture parameters on the permeability and permeability tensor were studied quantitatively.


## 2. Model Development

### 2.1. Physical Model

To quantitatively analyse the directional pressure-sensitive effect of fractures, it was necessary to clarify the stress condition of the fracture system. The conventional stress analysis of fracture deformation was mainly based on the constant-stress boundary condition. However, studies that focus on the stress-sensitive characteristics under constant-volume boundary condition are rare. This is because many studies ignore the pressure-arching effect in the reservoir [27], and it is difficult to realize the constant-volume-boundary stress sensitivity experiment by using the conventional core experiment method.Ignoring the pressure-arching effect for low permeability and tight reservoirs will exaggerate the influence of stress sensitivity effect on reservoir-development characteristics. This can lead to the incorrect formulation of development strategies for such reservoirs. When there is no stress-arching effect in a reservoir, and the pore-fluid pressure decreases, it is assumed that the vertical deformation of the reservoir (due to overlying stress) isΔZ. Then, the overlying strata in the strain-swept area will also produce a ΔZ deformation with the deformation of the reservoir (Figure 1(a)). When a stress-arch effect exists, it can prevent the overlying strata from further deformation. The stress-arch partially offsets the overlying stress, and may even decrease it to (near) zero. If it is assumed that the deformation of a reservoir with decreasing pore-fluid pressure is ΔZ1 (a stress-arch effect exists), then ΔZ1 will always be below ΔZ, and ΔZ1 may be approximately 0 (Figure 1(a)). Furthermore, the value of ΔZ1 is mainly related to the mechanical properties of the surrounding rock and reservoir (elastic modulus E, Poisson’s ratio V), the shape of the reservoir. In addition, it also reflects the stress on the reservoir. Regardless of the presence of a stress-arch effect, the effect of lateral strata on the reservoir is the same as for overlying strata (Figure 1(b)).Figure 1
Stress-arching effect on the reservoir boundary displacement.
(a)
Overlying strata action(b)
Lateral strata actionTherefore, when there is no effect of a stress arch or when a boundary displacement—which is caused by the total stress acting on the reservoir after the stress arch had been balanced—cannot be ignored, the outer boundary of the reservoir is a constant-stress boundary. In addition, when the boundary displacement—which is caused by the total stress acting on the reservoir after the stress arch is balanced—is very small (and can be approximately ignored), the boundary of the reservoir can be regarded as “constant-volume.”When the rock around the reservoir was very hard (mechanical properties are similar to a rigid body), this was equivalent to adding a layer of a “hard shell” to the reservoir. The hard overlying strata can more effectively bear the external load, and the stress arch is easier to produce. At this point in time, regardless of how much stress the surrounding strata was subjected to, the overlying strata will show no obvious deformation. In other words, the deformation of the overlying strataΔZ1≈0 and the total reservoir-volume remains unchanged (Figure 2).Figure 2
Schematic illustrating the reservoir boundaries.


### 2.2. The Number of Fractures

Let us take typical microelements which contain a group of parallel fractures as the research object (seen in the square ABCD in Figure3(a)), and let these microelements be collectively referred to as the representative element (RE), hereafter. As shown in Figure 3, the angle between the fracture direction and the positive direction of the x-axis is β, and the fractures are distributed evenly on the plane (the fracture distance is d). The width, length, and height of the RE are L1, L2, and h, respectively. The RE is surrounded by fixed boundaries. For the fractures above the red dashed line in Figure 3(a), the fractures intersect the side AB. For example, point E is the intersection of fracture 1 and side AB, and the distance from point A to point E is d1. However, for fractures under the red dashed line, the fractures intersect the side BC. For example, point F is the intersection of fracture 1′ and edge BC, and the distance from point C to point F is d2. Mathematically, d1 and d2 are the fracture distances along the x-axis and y-axis, respectively.Figure 3
Schematic diagram of fracture networks.
(a)
Single group of fractures(b)
Multiple groups of fracturesBased on Figure3(a), for fractured porous media with a single group of fractures, the number of fractures above the red line n1 is
(1)n1=L1d1=L1cosβd,wheredenotes a rounding operation. In addition, the number of fractures under the red line n2 is
(2)n2=L2d2=L2sinβd.With the assumptionL1=L2=L, the total number n of fractures in Figure 3(a) is
(3)n=n1+n2=Lsinβ+cosβd.Equation (3) reveals that n is the function of angle β, L, and d. Under the given L and d, n varies with angle β. For example, results (Table 1) suggest that, when fracture distance is set as 20 mm (d=d0=20mm) and L is set to 200 mm, the number of fractures varies from 10 to 14.14 as angle β increases from 0° to 45°. However, when β increases from 45° to 90°, the number of fractures decreases from 14.14 to 10. To retain the number of fractures, the fracture distance should be modified. Table 2 shows the quantitative relationship between the fracture distance and d0 (fracture distance for fractures with angle β=0°).Table 1
Number of fractures in the models with different fracture angles.
Fracture angle (°)Number of fractures (strip)Fracture angle (°)Number of fractures (strip)0n0=10.0050n50=14.095n5=10.8355n55=13.9310n10=11.5860n60=13.6615n15=12.2565n65=13.2920n20=12.8270n70=12.8225n25=13.2975n75=12.2530n30=13.6680n80=11.5835n35=13.9385n85=10.8340n40=14.0990n90=10.0045n45=14.14Table 2
Relationship between fracture distancedx and d0.
Fracture angle (°)Geometric relationship withd0Fracture angle (°)Geometric relationship withd00d050d50=2d0cos5°5d5=2d0cos40°55d55=2d0cos10°10d10=2d0cos35°60d60=1+3d0/215d15=6/2d065d65=2d0cos20°20d20=2d0cos25°70d70=2d0cos25°25d25=2d0cos20°75d75=6/2d030d30=1+3d0/280d80=2d0cos35°35d35=2d0cos10°85d85=2d0cos40°40d40=2d0cos5°90d045d45=2d0Based on equation (3), for fractured porous media with multiple groups of fractures (Figure 3(b)), the total number mt of fractures is
(4)mt=∑i=1mni=∑i=1mLisinβi+cosβidi.


### 2.3. Effective Stress of Fractured Porous Media

Based on the concept of effective stress in soil mechanics firstly proposed by Terzaghi [16], the effective stress at a given point can be determined as
(5)σei=σi−p,i=x,y,z,where σei denotes the effective stress in the ii=x,y,z direction, σi denotes the total stress (confining pressure) in the i direction, and p is the pore-fluid pressure.In 1961, according to the force balance on particles, Skempton [28] modified equation (5), which could be written as
(6)σe=σ−αp,where α is the effective-stress coefficient (Biot coefficient), which is
(7)α=1−CsC,where Cs/C denotes the ratio of the volume compressibility of rock particles to rock skeleton volume compressibility. Physically, α varies from 0 to 1, and it will be affected by various factors, such as the cementation way of the porous medium, rock lithology, porosity, and rock structure characterization. As reported in the literature [29], the ratio ranged from 0.04 to 0.5. In general, for unconsolidated or weak rocks, α is approximately 1.


### 2.4. Relationship between Displacement of Matrix System and Fracture System

Physically, fractured porous media are composed of the matrix system and the fracture system. Due to effective stress, both these two systems will deform. To simplify the model, the stress field is assumed to be isotropic and the deformation of the fractured porous media due to effective stress is instantaneous. Moreover, during the deformation, the displacement of the outer boundary of the porous media is zero. For a representative elementary (RE) of fractured porous media, we assume that both the initial length and initial width areL0, the initial thickness is h0, and the overall equivalent elastic modulus and Poisson’s ratio of RE are Et and νt, respectively.We posit the assumption that, even though the experimental model is affected by the effective stress, the total volume of the model remains unchanged (i.e., the displacement of the outer boundary is zero). Thus, we have [2, 30]
(8)Δux=Δumx+Δufx=L0EtΔσex−νtΔσey+Δσez=0,Δuy=Δumy+Δufy=L0EtΔσey−νtΔσex+Δσez=0,Δuz=Δumz=h0EtΔσez−νtΔσex+Δσey=0,where Δux, Δuy, and Δuz are the total displacement of the model along the x, y, and z directions, respectively; Δumx, Δumy, and Δumz are the total displacement of the matrix system along the x, y, and z directions, respectively; Δufx and Δufy are the total displacement of the fracture system along the x, y, and z directions, respectively; Δσex, Δσey, and Δσez are the effective-stress changes of the model along the x, y, and z directions, respectively. It should be noted that, the direction of compressive stress is positive and the direction of tensile stress is negative.By solving equation (8), we have
(9)Δσex=Δσey=Δσez.Based on effective-stress principle [28, 31, 32], the effective stress of the model in the three directions (x, y, and z) can be expressed as follows:
(10)Δσx−αΔp=Δσy−αΔp=Δσz−αΔp,where Δσx, Δσy, and Δσz are the total stress changes in the x, y, and z directions, respectively. α is the effective-stress coefficient of the model, and Δp is the change in pore-fluid pressure.Therefore,(11)Δσx=Δσy=Δσz.


### 2.5. Stress-Dependent Permeability of Fracture Networks

Substituting equation (A.6) in the appendix and equation (B.5) into equation (8), we obtain
(12)nd0sinβEm1−ν2Δσtmx−nd0sinβEmν+ν2Δσtmy+nd0sinβEm2ν2+ν−1αmΔpm+n⋅b0Δσx−αfΔpfEfsinβ=0,nd0cosβEm1−ν2Δσtmy−nd0cosβEmν+ν2Δσtmx+nd0cosβEm2ν2+ν−1αmΔpm+n⋅b0Δσy−αfΔpfEfcosβ=0.Here, the total stress change of both the matrix system and fracture system are equal, namely,Δσtmx=Δσx, Δσtmy=Δσy, Δσtmz=Δσz, and Δσx=Δσy=Δσz. In addition, we assume that the fluid pressure in the matrix system and the fracture system reaches equilibrium instantly, i.e., Δpm=Δpf=Δp and Δp=p−p0, where p denotes the current pore-fluid pressure.Solving equation (12), we obtain
(13)Δσx=Δσy=Emαfb0+d0Efαm1−2ν1+νEmb0+d0Ef1−2ν1+νΔp.Substituting equation (13) into equation (B.5) in the appendix, we can write
(14)Δufx=L0sinβ+cosβd0⋅b0Ef⋅d0Efαm−αf1−2ν1+νd0Ef1−2ν1+ν+b0EmΔp⋅sinβ,Δufy=L0sinβ+cosβd0⋅b0Ef⋅d0Efαm−αf1−2ν1+νd0Ef1−2ν1+ν+b0EmΔp⋅cosβ.Thus, the total change for the aperture in the fracture system is(15)Δun=Δb=Δufx2+Δufy2=±L0sinβ+cosβd0⋅b0Efd0Efαm−αf1−2ν1+νd0Ef1−2ν1+ν+b0EmΔp.WhenΔun is positive, the fracture aperture decreases, and when it is negative, the fracture aperture increases.Therefore, the total fracture aperture after deformation can be written as(16)b=b0+Δb=L0b0sinβ+cosβd0⋅1+1Ef⋅d0Efαf−αm1−2ν1+νd0Ef1−2ν1+ν+b0EmΔp.Assuming that the fracture permeability iskf0 when the pore-fluid pressure is p0 and kf when the pore-fluid pressure is p, then the relationship between these parameters and the fracture aperture is as follows [4, 33]:
(17)kfkf0=1+Δbb03.kf0 can be calculated as follows:
(18)kf0=L0sinβ+cosβb0d0L0⋅b0212=sinβ+cosβb0312d0.After substituting equation (18) and equation (15) into equation (17), the permeability of the fracture system in parallel fracture directions can be obtained:
(19)kf=b03sinβ+cosβ12d0⋅1+1Efd0Efαf−αm1−2ν1+νd0Ef1−2ν1+ν+b0EmΔp3.In fractured porous media,0<αm<1, 0<αf<1, and αf>αm. In addition, when the fractures in porous media do not develop, it is possible that αf=αm [3].


## 2.1. Physical Model

To quantitatively analyse the directional pressure-sensitive effect of fractures, it was necessary to clarify the stress condition of the fracture system. The conventional stress analysis of fracture deformation was mainly based on the constant-stress boundary condition. However, studies that focus on the stress-sensitive characteristics under constant-volume boundary condition are rare. This is because many studies ignore the pressure-arching effect in the reservoir [27], and it is difficult to realize the constant-volume-boundary stress sensitivity experiment by using the conventional core experiment method.Ignoring the pressure-arching effect for low permeability and tight reservoirs will exaggerate the influence of stress sensitivity effect on reservoir-development characteristics. This can lead to the incorrect formulation of development strategies for such reservoirs. When there is no stress-arching effect in a reservoir, and the pore-fluid pressure decreases, it is assumed that the vertical deformation of the reservoir (due to overlying stress) isΔZ. Then, the overlying strata in the strain-swept area will also produce a ΔZ deformation with the deformation of the reservoir (Figure 1(a)). When a stress-arch effect exists, it can prevent the overlying strata from further deformation. The stress-arch partially offsets the overlying stress, and may even decrease it to (near) zero. If it is assumed that the deformation of a reservoir with decreasing pore-fluid pressure is ΔZ1 (a stress-arch effect exists), then ΔZ1 will always be below ΔZ, and ΔZ1 may be approximately 0 (Figure 1(a)). Furthermore, the value of ΔZ1 is mainly related to the mechanical properties of the surrounding rock and reservoir (elastic modulus E, Poisson’s ratio V), the shape of the reservoir. In addition, it also reflects the stress on the reservoir. Regardless of the presence of a stress-arch effect, the effect of lateral strata on the reservoir is the same as for overlying strata (Figure 1(b)).Figure 1
Stress-arching effect on the reservoir boundary displacement.
(a)
Overlying strata action(b)
Lateral strata actionTherefore, when there is no effect of a stress arch or when a boundary displacement—which is caused by the total stress acting on the reservoir after the stress arch had been balanced—cannot be ignored, the outer boundary of the reservoir is a constant-stress boundary. In addition, when the boundary displacement—which is caused by the total stress acting on the reservoir after the stress arch is balanced—is very small (and can be approximately ignored), the boundary of the reservoir can be regarded as “constant-volume.”When the rock around the reservoir was very hard (mechanical properties are similar to a rigid body), this was equivalent to adding a layer of a “hard shell” to the reservoir. The hard overlying strata can more effectively bear the external load, and the stress arch is easier to produce. At this point in time, regardless of how much stress the surrounding strata was subjected to, the overlying strata will show no obvious deformation. In other words, the deformation of the overlying strataΔZ1≈0 and the total reservoir-volume remains unchanged (Figure 2).Figure 2
Schematic illustrating the reservoir boundaries.


## 2.2. The Number of Fractures

Let us take typical microelements which contain a group of parallel fractures as the research object (seen in the square ABCD in Figure3(a)), and let these microelements be collectively referred to as the representative element (RE), hereafter. As shown in Figure 3, the angle between the fracture direction and the positive direction of the x-axis is β, and the fractures are distributed evenly on the plane (the fracture distance is d). The width, length, and height of the RE are L1, L2, and h, respectively. The RE is surrounded by fixed boundaries. For the fractures above the red dashed line in Figure 3(a), the fractures intersect the side AB. For example, point E is the intersection of fracture 1 and side AB, and the distance from point A to point E is d1. However, for fractures under the red dashed line, the fractures intersect the side BC. For example, point F is the intersection of fracture 1′ and edge BC, and the distance from point C to point F is d2. Mathematically, d1 and d2 are the fracture distances along the x-axis and y-axis, respectively.Figure 3
Schematic diagram of fracture networks.
(a)
Single group of fractures(b)
Multiple groups of fracturesBased on Figure3(a), for fractured porous media with a single group of fractures, the number of fractures above the red line n1 is
(1)n1=L1d1=L1cosβd,wheredenotes a rounding operation. In addition, the number of fractures under the red line n2 is
(2)n2=L2d2=L2sinβd.With the assumptionL1=L2=L, the total number n of fractures in Figure 3(a) is
(3)n=n1+n2=Lsinβ+cosβd.Equation (3) reveals that n is the function of angle β, L, and d. Under the given L and d, n varies with angle β. For example, results (Table 1) suggest that, when fracture distance is set as 20 mm (d=d0=20mm) and L is set to 200 mm, the number of fractures varies from 10 to 14.14 as angle β increases from 0° to 45°. However, when β increases from 45° to 90°, the number of fractures decreases from 14.14 to 10. To retain the number of fractures, the fracture distance should be modified. Table 2 shows the quantitative relationship between the fracture distance and d0 (fracture distance for fractures with angle β=0°).Table 1
Number of fractures in the models with different fracture angles.
Fracture angle (°)Number of fractures (strip)Fracture angle (°)Number of fractures (strip)0n0=10.0050n50=14.095n5=10.8355n55=13.9310n10=11.5860n60=13.6615n15=12.2565n65=13.2920n20=12.8270n70=12.8225n25=13.2975n75=12.2530n30=13.6680n80=11.5835n35=13.9385n85=10.8340n40=14.0990n90=10.0045n45=14.14Table 2
Relationship between fracture distancedx and d0.
Fracture angle (°)Geometric relationship withd0Fracture angle (°)Geometric relationship withd00d050d50=2d0cos5°5d5=2d0cos40°55d55=2d0cos10°10d10=2d0cos35°60d60=1+3d0/215d15=6/2d065d65=2d0cos20°20d20=2d0cos25°70d70=2d0cos25°25d25=2d0cos20°75d75=6/2d030d30=1+3d0/280d80=2d0cos35°35d35=2d0cos10°85d85=2d0cos40°40d40=2d0cos5°90d045d45=2d0Based on equation (3), for fractured porous media with multiple groups of fractures (Figure 3(b)), the total number mt of fractures is
(4)mt=∑i=1mni=∑i=1mLisinβi+cosβidi.


## 2.3. Effective Stress of Fractured Porous Media

Based on the concept of effective stress in soil mechanics firstly proposed by Terzaghi [16], the effective stress at a given point can be determined as
(5)σei=σi−p,i=x,y,z,where σei denotes the effective stress in the ii=x,y,z direction, σi denotes the total stress (confining pressure) in the i direction, and p is the pore-fluid pressure.In 1961, according to the force balance on particles, Skempton [28] modified equation (5), which could be written as
(6)σe=σ−αp,where α is the effective-stress coefficient (Biot coefficient), which is
(7)α=1−CsC,where Cs/C denotes the ratio of the volume compressibility of rock particles to rock skeleton volume compressibility. Physically, α varies from 0 to 1, and it will be affected by various factors, such as the cementation way of the porous medium, rock lithology, porosity, and rock structure characterization. As reported in the literature [29], the ratio ranged from 0.04 to 0.5. In general, for unconsolidated or weak rocks, α is approximately 1.


## 2.4. Relationship between Displacement of Matrix System and Fracture System

Physically, fractured porous media are composed of the matrix system and the fracture system. Due to effective stress, both these two systems will deform. To simplify the model, the stress field is assumed to be isotropic and the deformation of the fractured porous media due to effective stress is instantaneous. Moreover, during the deformation, the displacement of the outer boundary of the porous media is zero. For a representative elementary (RE) of fractured porous media, we assume that both the initial length and initial width areL0, the initial thickness is h0, and the overall equivalent elastic modulus and Poisson’s ratio of RE are Et and νt, respectively.We posit the assumption that, even though the experimental model is affected by the effective stress, the total volume of the model remains unchanged (i.e., the displacement of the outer boundary is zero). Thus, we have [2, 30]
(8)Δux=Δumx+Δufx=L0EtΔσex−νtΔσey+Δσez=0,Δuy=Δumy+Δufy=L0EtΔσey−νtΔσex+Δσez=0,Δuz=Δumz=h0EtΔσez−νtΔσex+Δσey=0,where Δux, Δuy, and Δuz are the total displacement of the model along the x, y, and z directions, respectively; Δumx, Δumy, and Δumz are the total displacement of the matrix system along the x, y, and z directions, respectively; Δufx and Δufy are the total displacement of the fracture system along the x, y, and z directions, respectively; Δσex, Δσey, and Δσez are the effective-stress changes of the model along the x, y, and z directions, respectively. It should be noted that, the direction of compressive stress is positive and the direction of tensile stress is negative.By solving equation (8), we have
(9)Δσex=Δσey=Δσez.Based on effective-stress principle [28, 31, 32], the effective stress of the model in the three directions (x, y, and z) can be expressed as follows:
(10)Δσx−αΔp=Δσy−αΔp=Δσz−αΔp,where Δσx, Δσy, and Δσz are the total stress changes in the x, y, and z directions, respectively. α is the effective-stress coefficient of the model, and Δp is the change in pore-fluid pressure.Therefore,(11)Δσx=Δσy=Δσz.


## 2.5. Stress-Dependent Permeability of Fracture Networks

Substituting equation (A.6) in the appendix and equation (B.5) into equation (8), we obtain
(12)nd0sinβEm1−ν2Δσtmx−nd0sinβEmν+ν2Δσtmy+nd0sinβEm2ν2+ν−1αmΔpm+n⋅b0Δσx−αfΔpfEfsinβ=0,nd0cosβEm1−ν2Δσtmy−nd0cosβEmν+ν2Δσtmx+nd0cosβEm2ν2+ν−1αmΔpm+n⋅b0Δσy−αfΔpfEfcosβ=0.Here, the total stress change of both the matrix system and fracture system are equal, namely,Δσtmx=Δσx, Δσtmy=Δσy, Δσtmz=Δσz, and Δσx=Δσy=Δσz. In addition, we assume that the fluid pressure in the matrix system and the fracture system reaches equilibrium instantly, i.e., Δpm=Δpf=Δp and Δp=p−p0, where p denotes the current pore-fluid pressure.Solving equation (12), we obtain
(13)Δσx=Δσy=Emαfb0+d0Efαm1−2ν1+νEmb0+d0Ef1−2ν1+νΔp.Substituting equation (13) into equation (B.5) in the appendix, we can write
(14)Δufx=L0sinβ+cosβd0⋅b0Ef⋅d0Efαm−αf1−2ν1+νd0Ef1−2ν1+ν+b0EmΔp⋅sinβ,Δufy=L0sinβ+cosβd0⋅b0Ef⋅d0Efαm−αf1−2ν1+νd0Ef1−2ν1+ν+b0EmΔp⋅cosβ.Thus, the total change for the aperture in the fracture system is(15)Δun=Δb=Δufx2+Δufy2=±L0sinβ+cosβd0⋅b0Efd0Efαm−αf1−2ν1+νd0Ef1−2ν1+ν+b0EmΔp.WhenΔun is positive, the fracture aperture decreases, and when it is negative, the fracture aperture increases.Therefore, the total fracture aperture after deformation can be written as(16)b=b0+Δb=L0b0sinβ+cosβd0⋅1+1Ef⋅d0Efαf−αm1−2ν1+νd0Ef1−2ν1+ν+b0EmΔp.Assuming that the fracture permeability iskf0 when the pore-fluid pressure is p0 and kf when the pore-fluid pressure is p, then the relationship between these parameters and the fracture aperture is as follows [4, 33]:
(17)kfkf0=1+Δbb03.kf0 can be calculated as follows:
(18)kf0=L0sinβ+cosβb0d0L0⋅b0212=sinβ+cosβb0312d0.After substituting equation (18) and equation (15) into equation (17), the permeability of the fracture system in parallel fracture directions can be obtained:
(19)kf=b03sinβ+cosβ12d0⋅1+1Efd0Efαf−αm1−2ν1+νd0Ef1−2ν1+ν+b0EmΔp3.In fractured porous media,0<αm<1, 0<αf<1, and αf>αm. In addition, when the fractures in porous media do not develop, it is possible that αf=αm [3].


## 3. Results and Discussions

### 3.1. Model Validation

#### 3.1.1. Experimental Observations

To better observe the dynamic change process of the fracture aperture when the fluid pressure changes, a transparent plexiglass tank was used as the container in the experimental setup (Figure4). The measurement setup mainly consisted of three parts (Figure 5): fluid-injection system (constant-pressure and constant-flow pump, intermediate container), fluid-displacement system (six-way valve, fracture-seepage model, and back-pressure valve), and measurement system (pressure gauge and electronic balance).Figure 4
Plexiglass container.Figure 5
Schematic of the measurement setup.As shown in Figure6, nine angle-models (β=0°, 10°, 20°, 30°, 45°, 60°, 70°, 80°, and 90°, respectively) were used. The measurement process mainly includes the following two steps: (1) Simultaneous adjustments of the pressure at the inlet and outlet of the model are performed to keep the pressure gradient constant. Changes are made to the average pore-fluid pressure, and the pressure at both ends of the model and the outlet flow at the outlet end for different average pore-fluid pressures are recorded. (2) Simultaneous adjustments of the inlet and outlet pressures of the model are performed to keep the average pore-fluid pressure constant. Changes are made to the pressure gradient, and the pressure at both ends of the model and the flow at the outlet end for different pressure gradient conditions are recorded.Figure 6
Schematic diagrams of different fracture direction model designs.When the pressures at both ends of the entrance and exit are 0 MPa, the average pressure of the model is 0 MPa, the pressure gradient is 0 MPa/m, and the fracture is closed. When pressures at both inlet and outlet increases from 0 MPa (Figure7(a)) to 0.016 MPa (Figure 7(b)), the fracture is open. As a result, the inlet and outlet pressures continue to increase from 0.016 MPa to 0.025 MPa, and the fracture aperture further increases significantly (Figure 7(c)). Furthermore, it can be seen that the fracture aperture in the model increases with increasing pore-fluid pressure.Figure 7
Fracture aperture as a function of average pore-fluid pressure.
(a)
0 MPa(b)
0.016 MPa(c)
0.025 MPaAt the initial state of the model, the pressure gradient was 0 MPa/m, the average pore-fluid pressure was 0 MPa, and the fracture was closed (Figure8(a)). During injection and production, along the direction of the pressure gradient, the fracture changed from ||-shaped to V-shaped. In addition, the distribution of fracture width was uneven. Furthermore, the fracture width near the inlet end was large, but it was small near the outlet (Figure 8(b)). This was the case because the pressure at the inlet was high, which means the fracture width was large. On the other hand, the pressure at the outlet was small, which caused the fracture width to be small. Figure 8(c) is a magnification of the fracture at the top of the model in Figure 8(b). There are clear differences in the fracture aperture between the inlet end and the outlet. When the outlet pressure remains unchanged and the pressure gradient is increased by increasing the inlet pressure, the average fracture aperture increased. However, when the inlet pressure remained unchanged, and the pressure gradient is increased by decreasing the outlet pressure, the average fracture aperture decreased.Figure 8
Average fracture aperture as a function of the pressure gradient.
(a)(b)(c)The experimental results indicate that the deformation characteristics of fractures in the model are not the same before and after water injection, while the average aperture of fractures was affected by the uneven pressure at the inlet and outlet.


#### 3.1.2. Consistency between Theoretical and Experimental Results

The equation for motion of the two-dimensional physical simulation model is [34] as follows:
(20)V⟶=K⟶⋅∂p∂x∂p∂y=kxx∂p∂x+kxy∂p∂ykyx∂p∂x+kyy∂p∂y=vxvy.Because each experimental test point was in a steady flow state, the governing equation is as follows:(21)kxx∂2p∂x2+kxy∂2p∂x∂y+kyx∂2p∂y∂x+kyy∂2p∂y2=0,and the boundary conditions are as follows:
(22)p=p1,x=0,y∈0,L,p=p2,x=L,y∈0,L,vy=0,y=L,x∈0,L,vy=0,y=0,x∈0,L.Here,kxx, kxy, kyx, and kyy denote the components of the permeability tensor in each direction, where kxy=kyx; ∂p/∂x and ∂p/∂y are the pressure gradients in the x and y directions, respectively; V⟶ is the velocity vector; vx and vy are the seepage velocities in the x and y directions, respectively; p1 and p2 denote the pressure at the inlet and outlet, respectively; the length and width of the model are defined by L.The governing equation shown in equation (21) is difficult to use to directly calculate the different components of the permeability tensor. Therefore, with the help of the Schwarz-Christoffel transformation [35] (Figure 9), the flow calculation equation could be performed.Figure 9
The Schwarz-Christoffel map that maps the area inside the polygon to a rectangle.If fracture permeabilitykf is known, flow rate can be calculated using the following formula:
(23)Q=kf+ky⋅ky1/2⋅Δp⋅hμ⋅KmK′m.Here,Q denotes the flow at the model exit; ky is permeability in y-direction; μ is fluid viscosity; h is model thickness; Km is the first type of complete elliptic integral with parameter m, where K′m is the first type of complete elliptic integral with parameter 1−m [36]. Furthermore, Km and K′m correspond to a specific value of g, and the value of Km and K′m can be found in the complete elliptic integral table [37]. Alternatively, it can be calculated via numerical modeling. The solution procedure for equation (23) is shown in Figure 10.Figure 10
Fracture permeability semianalytical model solver calculation flow chart.We considered the experimental results of the 0° fracture model as a sample. There were ten fractures in the model, where the angle between all fractures andx-axis was 0°, and the fracture distance was 2 cm. The basic parameters of the model are shown in Table 3. It can be seen from Figure 11 that the experimental test data are in good agreement with the theoretical calculation results (equation (23)). This suggests that the fracture permeability calculation model used in this paper is sufficiently accurate.Table 3
Summary of the parameters used in the model.
Physical parametersValueUnitPhysical parametersValueUnitInitial matrix porosity,φm05%Elastic modulus of fracture system,Ef0.5MPaInitial matrix permeability,km00.4μm2Fluid viscosity,μ1mPa·sInitial fracture permeability,kf00.92μm2Model length, model width,L200mmElastic modulus of matrix system,Em5.112MPaModel thickness,h30mmInitial fracture aperture,b00.07mmEffective-stress coefficient of matrix system,αm0.65—Fracture distance,d020mmEffective-stress coefficient of fracture system,αf1—Figure 11
Comparison between experimental and theoretical outcomes.
(a)
Flow rate(b)
Fracture permeability


### 3.2. Influencing Factors of Fracture Pressure Sensitivity Characteristics

The parameters, which affect the variation of fracture aperture, are mainly the initial fracture aperture, pore-fluid pressure, fracture distance, elastic modulus, effective-stress coefficient, and the fracture direction. We will use the new pressure-sensitivity equation (equation (19)) to study fracture aperture and permeability change with these parameters. The most important parameters that were used for the model calculation are shown in Table 3.


#### 3.2.1. Initial Fracture Apertureb0

The pore-fluid pressure remained constant. The larger the initial aperture of the fracture was, the larger was the aperture and corresponding fracture permeability following deformation (Figure12).Figure 12
Fracture aperture and fracture permeability vs. initial fracture aperture.


#### 3.2.2. Pore-Fluid Pressurep

As shown in Figure13, the initial fracture apertures were identical. Both, the fracture aperture and fracture permeability increased with increasing of fluid pressure.Figure 13
Fracture aperture and fracture permeability vs. pore-fluid pressure.


#### 3.2.3. Initial Fracture Distanced0

The larger the fracture distance, the smaller the fracture density and the smaller the number of fractures in the same size unit. Therefore, when fluid pressure remained constant, the larger the fracture distance, the smaller the total apertures and the smaller the corresponding fracture permeability (Figure14).Figure 14
Fracture aperture and fracture permeability vs. initial fracture space.


#### 3.2.4. Elastic Modulus

The elastic modulus of the matrix and fracture systems are related to each other. For the same seepage medium, it is necessary to discuss the change law for fracture permeability, when the two variables change simultaneously. As shown in Figure15, when the changes in effective stress are identical, the greater is the elastic modulus of the matrix system, and the smaller is the strain for the matrix system. Because the displacement values for the fracture system and the matrix system are equal in magnitude and opposite in direction, the deformation of the fracture system is smaller. Therefore, both fracture aperture and permeability decrease with increasing elastic modulus of the matrix system. When the elastic modulus of the fracture system increased, the change law did the same.Figure 15
Fracture aperture and fracture permeability vs. fracture modulus.


#### 3.2.5. Effective-Stress Coefficient

There is also a specific relationship between the effective-stress coefficient of the matrix and the fracture system. However, there is no situation where one changes while the other remains unchanged. Therefore, for the same RE, it is necessary to discuss the change law of fracture permeability, when they change simultaneously. As shown in Figure16, the larger the effective-stress coefficient was, the smaller were the effective stress and the deformation of the matrix and fracture under the same total stress change and pore-fluid change. The larger the effective-stress coefficient was—following after the fracture has deformed—the smaller were the aperture and permeability.Figure 16
Fracture aperture and fracture permeability vs. effective-stress coefficient.


#### 3.2.6. Fracture Angleβ

As the fracture angleβ increased, the total number of fractures in RE first increased and then decreased. Therefore, both the total fracture aperture and fracture permeability in isotropic media increase first and then decrease.Whenβ was 45°, the fracture permeability reached the maximum. When β was 0° and 90°, the fracture permeability was equal and reached the minimum (Figure 17).Figure 17
Fracture aperture and fracture permeability with fracture direction.


### 3.3. Influencing Factors for the Fracture Permeability Tensor

We introduced two rectangular coordinate systems (x′y′, xy) along the direction of fracture development and the direction of macroscopic pressure gradient, respectively. Assuming matrix permeability was km, the permeability tensor of RE in a fracture coordinate system is as follows:
(24)K′⟶=kf+km00km.After converting equation (24) to the macroscopic pressure-gradient coordinate system xy via coordinate transformation, the associated permeability tensor becomes [19, 38]
(25)K⟶=kxxkxykyxkyy=kf+2km2+kf2cos2βkf2sin2βkf2sin2βkf+2km2−kf2cos2β.Becausekm≪kf, equation (25) can be simplified to:
(26)K⟶=kfcos2βsinβcosβsinβcosβsin2β.The angle between the direction of the main permeability and thex-axis is as follows:
(27)β=12arctan2kxykxx−kxy+π2,β≤−π4,β=12arctan2kxykxx−kxy,−π4≤β≤π4,β=12arctan2kxykxx−kxy−π2,β≥π4.


#### 3.3.1. Permeability-Tensor Componentskxx, kxy, and kyy Vary with Fracture Angle

The permeability-tensor componentkxx decreases with increasing fracture angle, kyy increases with increasing fracture angle, and kxy increases first but then decreases, with increasing fracture angle. When the fracture angle was 45°, the component kxy can reach its maximum value. When the fracture angle was 0° and 90°, the component kxy was 0 (Figure 18). It shows that as the fracture angle increased, the flow capacity in the x direction weakened gradually, while the fluid-flow capacity in the y direction gradually increased.Figure 18
Variation of the permeability tensor with fracture direction.


#### 3.3.2. Variation Characteristic of Main Value Direction of Permeability

Within an isotropic seepage medium without fractures, the permeability was the same everywhere in the medium. According to equation (25), the permeability tensor ellipse of fractured media can be obtained (Figure 19). When a group of fractures develops in the media, the overall permeability of the media can show anisotropic characteristics. Moreover, the maximum main value direction for permeability was always along the fracture-development direction, while the minimum principal value direction was always perpendicular to the fracture-development direction.Figure 19
Anisotropic permeability tensor vs. fracture angle change.


#### 3.3.3. Influence of Fracture Deformation on Permeability Tensor

Figure20(a) shows the variation of a permeability tensor vs. pore-fluid pressure in a seepage medium with a group of 0° fracture development. Figure 20(b) shows the variation of a permeability tensor with pore-fluid pressure in a seepage medium for a group of 45° fracture development. With increasing fluid pressure, the maximum principal value of medium permeability (length of elliptical semimajor axis) increases. However, the main value direction remains unchanged, and vice versa. This suggests that, given a seepage medium with a set of fracture development and isotropic matrix, the fracture deformation only alters the main value of the permeability.Figure 20
Anisotropic permeability tensor as a function of fluid pressure.
(a)
A set of 0° fractures(b)
A set of 45° fracturesUsing the two groups of fracture distributions shown in Figure21 as an example, when two groups of unsymmetrical fractures develop in RE (Figure 21(a)), the principal value of permeability increased with increasing fluid pressure. In addition, the direction of the permeability principal value gradually approached the group of fractures with strong pressure sensitivity (Figure 22(a)). Moreover, when two groups of fractures follow an orthogonal distribution (Figure 21(b)), or the sum of two groups of fracture angles was 90°, only the principal value changed (Figure 22(b)).Figure 21
Schematic diagram for two groups of fracture models.
(a)
0°~60°(b)
0°~90°Figure 22
Variation of main values of fracture permeability for two groups at different angles.
(a)
0°~60°(b)
0°~90°


### 3.4. Mechanism and Theoretical Explanation of Tensor Permeability Variation

According to the analysis results in Sections3.1, 3.2, and 3.3, under pressure sensitivity, not only will the principal value of the full tensor permeability of fractured anisotropic media change but the direction of principal permeability will also rotate. With increasing injection-production pressure difference or decreasing average pressure, the principal value direction will be inclined to the development direction of weakly sensitive fractures, and vice versa. The mechanism of this change is that several groups of fractures in the anisotropic strain medium produce nonlinear deformation on different degrees under pressure-sensitive conditions, which destroys the ratio relationship of the original permeability tensor.From the perspective of tensor theory, the permeability field is a symmetric tensor of second order. The change of a permeability tensor field caused by pressure sensitivity will inevitably lead to a “transformation effect” in the coordinate system of the permeability tensor. The transformation effect of the coordinate system has and only includes two aspects: (1) the expansion of the coordinate size, that is, the change of the principal value of permeability, and (2) the rotation of the coordinate system, that is, the change of the principal direction of permeability.


## 3.1. Model Validation

### 3.1.1. Experimental Observations

To better observe the dynamic change process of the fracture aperture when the fluid pressure changes, a transparent plexiglass tank was used as the container in the experimental setup (Figure4). The measurement setup mainly consisted of three parts (Figure 5): fluid-injection system (constant-pressure and constant-flow pump, intermediate container), fluid-displacement system (six-way valve, fracture-seepage model, and back-pressure valve), and measurement system (pressure gauge and electronic balance).Figure 4
Plexiglass container.Figure 5
Schematic of the measurement setup.As shown in Figure6, nine angle-models (β=0°, 10°, 20°, 30°, 45°, 60°, 70°, 80°, and 90°, respectively) were used. The measurement process mainly includes the following two steps: (1) Simultaneous adjustments of the pressure at the inlet and outlet of the model are performed to keep the pressure gradient constant. Changes are made to the average pore-fluid pressure, and the pressure at both ends of the model and the outlet flow at the outlet end for different average pore-fluid pressures are recorded. (2) Simultaneous adjustments of the inlet and outlet pressures of the model are performed to keep the average pore-fluid pressure constant. Changes are made to the pressure gradient, and the pressure at both ends of the model and the flow at the outlet end for different pressure gradient conditions are recorded.Figure 6
Schematic diagrams of different fracture direction model designs.When the pressures at both ends of the entrance and exit are 0 MPa, the average pressure of the model is 0 MPa, the pressure gradient is 0 MPa/m, and the fracture is closed. When pressures at both inlet and outlet increases from 0 MPa (Figure7(a)) to 0.016 MPa (Figure 7(b)), the fracture is open. As a result, the inlet and outlet pressures continue to increase from 0.016 MPa to 0.025 MPa, and the fracture aperture further increases significantly (Figure 7(c)). Furthermore, it can be seen that the fracture aperture in the model increases with increasing pore-fluid pressure.Figure 7
Fracture aperture as a function of average pore-fluid pressure.
(a)
0 MPa(b)
0.016 MPa(c)
0.025 MPaAt the initial state of the model, the pressure gradient was 0 MPa/m, the average pore-fluid pressure was 0 MPa, and the fracture was closed (Figure8(a)). During injection and production, along the direction of the pressure gradient, the fracture changed from ||-shaped to V-shaped. In addition, the distribution of fracture width was uneven. Furthermore, the fracture width near the inlet end was large, but it was small near the outlet (Figure 8(b)). This was the case because the pressure at the inlet was high, which means the fracture width was large. On the other hand, the pressure at the outlet was small, which caused the fracture width to be small. Figure 8(c) is a magnification of the fracture at the top of the model in Figure 8(b). There are clear differences in the fracture aperture between the inlet end and the outlet. When the outlet pressure remains unchanged and the pressure gradient is increased by increasing the inlet pressure, the average fracture aperture increased. However, when the inlet pressure remained unchanged, and the pressure gradient is increased by decreasing the outlet pressure, the average fracture aperture decreased.Figure 8
Average fracture aperture as a function of the pressure gradient.
(a)(b)(c)The experimental results indicate that the deformation characteristics of fractures in the model are not the same before and after water injection, while the average aperture of fractures was affected by the uneven pressure at the inlet and outlet.


### 3.1.2. Consistency between Theoretical and Experimental Results

The equation for motion of the two-dimensional physical simulation model is [34] as follows:
(20)V⟶=K⟶⋅∂p∂x∂p∂y=kxx∂p∂x+kxy∂p∂ykyx∂p∂x+kyy∂p∂y=vxvy.Because each experimental test point was in a steady flow state, the governing equation is as follows:(21)kxx∂2p∂x2+kxy∂2p∂x∂y+kyx∂2p∂y∂x+kyy∂2p∂y2=0,and the boundary conditions are as follows:
(22)p=p1,x=0,y∈0,L,p=p2,x=L,y∈0,L,vy=0,y=L,x∈0,L,vy=0,y=0,x∈0,L.Here,kxx, kxy, kyx, and kyy denote the components of the permeability tensor in each direction, where kxy=kyx; ∂p/∂x and ∂p/∂y are the pressure gradients in the x and y directions, respectively; V⟶ is the velocity vector; vx and vy are the seepage velocities in the x and y directions, respectively; p1 and p2 denote the pressure at the inlet and outlet, respectively; the length and width of the model are defined by L.The governing equation shown in equation (21) is difficult to use to directly calculate the different components of the permeability tensor. Therefore, with the help of the Schwarz-Christoffel transformation [35] (Figure 9), the flow calculation equation could be performed.Figure 9
The Schwarz-Christoffel map that maps the area inside the polygon to a rectangle.If fracture permeabilitykf is known, flow rate can be calculated using the following formula:
(23)Q=kf+ky⋅ky1/2⋅Δp⋅hμ⋅KmK′m.Here,Q denotes the flow at the model exit; ky is permeability in y-direction; μ is fluid viscosity; h is model thickness; Km is the first type of complete elliptic integral with parameter m, where K′m is the first type of complete elliptic integral with parameter 1−m [36]. Furthermore, Km and K′m correspond to a specific value of g, and the value of Km and K′m can be found in the complete elliptic integral table [37]. Alternatively, it can be calculated via numerical modeling. The solution procedure for equation (23) is shown in Figure 10.Figure 10
Fracture permeability semianalytical model solver calculation flow chart.We considered the experimental results of the 0° fracture model as a sample. There were ten fractures in the model, where the angle between all fractures andx-axis was 0°, and the fracture distance was 2 cm. The basic parameters of the model are shown in Table 3. It can be seen from Figure 11 that the experimental test data are in good agreement with the theoretical calculation results (equation (23)). This suggests that the fracture permeability calculation model used in this paper is sufficiently accurate.Table 3
Summary of the parameters used in the model.
Physical parametersValueUnitPhysical parametersValueUnitInitial matrix porosity,φm05%Elastic modulus of fracture system,Ef0.5MPaInitial matrix permeability,km00.4μm2Fluid viscosity,μ1mPa·sInitial fracture permeability,kf00.92μm2Model length, model width,L200mmElastic modulus of matrix system,Em5.112MPaModel thickness,h30mmInitial fracture aperture,b00.07mmEffective-stress coefficient of matrix system,αm0.65—Fracture distance,d020mmEffective-stress coefficient of fracture system,αf1—Figure 11
Comparison between experimental and theoretical outcomes.
(a)
Flow rate(b)
Fracture permeability


## 3.1.1. Experimental Observations

To better observe the dynamic change process of the fracture aperture when the fluid pressure changes, a transparent plexiglass tank was used as the container in the experimental setup (Figure4). The measurement setup mainly consisted of three parts (Figure 5): fluid-injection system (constant-pressure and constant-flow pump, intermediate container), fluid-displacement system (six-way valve, fracture-seepage model, and back-pressure valve), and measurement system (pressure gauge and electronic balance).Figure 4
Plexiglass container.Figure 5
Schematic of the measurement setup.As shown in Figure6, nine angle-models (β=0°, 10°, 20°, 30°, 45°, 60°, 70°, 80°, and 90°, respectively) were used. The measurement process mainly includes the following two steps: (1) Simultaneous adjustments of the pressure at the inlet and outlet of the model are performed to keep the pressure gradient constant. Changes are made to the average pore-fluid pressure, and the pressure at both ends of the model and the outlet flow at the outlet end for different average pore-fluid pressures are recorded. (2) Simultaneous adjustments of the inlet and outlet pressures of the model are performed to keep the average pore-fluid pressure constant. Changes are made to the pressure gradient, and the pressure at both ends of the model and the flow at the outlet end for different pressure gradient conditions are recorded.Figure 6
Schematic diagrams of different fracture direction model designs.When the pressures at both ends of the entrance and exit are 0 MPa, the average pressure of the model is 0 MPa, the pressure gradient is 0 MPa/m, and the fracture is closed. When pressures at both inlet and outlet increases from 0 MPa (Figure7(a)) to 0.016 MPa (Figure 7(b)), the fracture is open. As a result, the inlet and outlet pressures continue to increase from 0.016 MPa to 0.025 MPa, and the fracture aperture further increases significantly (Figure 7(c)). Furthermore, it can be seen that the fracture aperture in the model increases with increasing pore-fluid pressure.Figure 7
Fracture aperture as a function of average pore-fluid pressure.
(a)
0 MPa(b)
0.016 MPa(c)
0.025 MPaAt the initial state of the model, the pressure gradient was 0 MPa/m, the average pore-fluid pressure was 0 MPa, and the fracture was closed (Figure8(a)). During injection and production, along the direction of the pressure gradient, the fracture changed from ||-shaped to V-shaped. In addition, the distribution of fracture width was uneven. Furthermore, the fracture width near the inlet end was large, but it was small near the outlet (Figure 8(b)). This was the case because the pressure at the inlet was high, which means the fracture width was large. On the other hand, the pressure at the outlet was small, which caused the fracture width to be small. Figure 8(c) is a magnification of the fracture at the top of the model in Figure 8(b). There are clear differences in the fracture aperture between the inlet end and the outlet. When the outlet pressure remains unchanged and the pressure gradient is increased by increasing the inlet pressure, the average fracture aperture increased. However, when the inlet pressure remained unchanged, and the pressure gradient is increased by decreasing the outlet pressure, the average fracture aperture decreased.Figure 8
Average fracture aperture as a function of the pressure gradient.
(a)(b)(c)The experimental results indicate that the deformation characteristics of fractures in the model are not the same before and after water injection, while the average aperture of fractures was affected by the uneven pressure at the inlet and outlet.


## 3.1.2. Consistency between Theoretical and Experimental Results

The equation for motion of the two-dimensional physical simulation model is [34] as follows:
(20)V⟶=K⟶⋅∂p∂x∂p∂y=kxx∂p∂x+kxy∂p∂ykyx∂p∂x+kyy∂p∂y=vxvy.Because each experimental test point was in a steady flow state, the governing equation is as follows:(21)kxx∂2p∂x2+kxy∂2p∂x∂y+kyx∂2p∂y∂x+kyy∂2p∂y2=0,and the boundary conditions are as follows:
(22)p=p1,x=0,y∈0,L,p=p2,x=L,y∈0,L,vy=0,y=L,x∈0,L,vy=0,y=0,x∈0,L.Here,kxx, kxy, kyx, and kyy denote the components of the permeability tensor in each direction, where kxy=kyx; ∂p/∂x and ∂p/∂y are the pressure gradients in the x and y directions, respectively; V⟶ is the velocity vector; vx and vy are the seepage velocities in the x and y directions, respectively; p1 and p2 denote the pressure at the inlet and outlet, respectively; the length and width of the model are defined by L.The governing equation shown in equation (21) is difficult to use to directly calculate the different components of the permeability tensor. Therefore, with the help of the Schwarz-Christoffel transformation [35] (Figure 9), the flow calculation equation could be performed.Figure 9
The Schwarz-Christoffel map that maps the area inside the polygon to a rectangle.If fracture permeabilitykf is known, flow rate can be calculated using the following formula:
(23)Q=kf+ky⋅ky1/2⋅Δp⋅hμ⋅KmK′m.Here,Q denotes the flow at the model exit; ky is permeability in y-direction; μ is fluid viscosity; h is model thickness; Km is the first type of complete elliptic integral with parameter m, where K′m is the first type of complete elliptic integral with parameter 1−m [36]. Furthermore, Km and K′m correspond to a specific value of g, and the value of Km and K′m can be found in the complete elliptic integral table [37]. Alternatively, it can be calculated via numerical modeling. The solution procedure for equation (23) is shown in Figure 10.Figure 10
Fracture permeability semianalytical model solver calculation flow chart.We considered the experimental results of the 0° fracture model as a sample. There were ten fractures in the model, where the angle between all fractures andx-axis was 0°, and the fracture distance was 2 cm. The basic parameters of the model are shown in Table 3. It can be seen from Figure 11 that the experimental test data are in good agreement with the theoretical calculation results (equation (23)). This suggests that the fracture permeability calculation model used in this paper is sufficiently accurate.Table 3
Summary of the parameters used in the model.
Physical parametersValueUnitPhysical parametersValueUnitInitial matrix porosity,φm05%Elastic modulus of fracture system,Ef0.5MPaInitial matrix permeability,km00.4μm2Fluid viscosity,μ1mPa·sInitial fracture permeability,kf00.92μm2Model length, model width,L200mmElastic modulus of matrix system,Em5.112MPaModel thickness,h30mmInitial fracture aperture,b00.07mmEffective-stress coefficient of matrix system,αm0.65—Fracture distance,d020mmEffective-stress coefficient of fracture system,αf1—Figure 11
Comparison between experimental and theoretical outcomes.
(a)
Flow rate(b)
Fracture permeability


## 3.2. Influencing Factors of Fracture Pressure Sensitivity Characteristics

The parameters, which affect the variation of fracture aperture, are mainly the initial fracture aperture, pore-fluid pressure, fracture distance, elastic modulus, effective-stress coefficient, and the fracture direction. We will use the new pressure-sensitivity equation (equation (19)) to study fracture aperture and permeability change with these parameters. The most important parameters that were used for the model calculation are shown in Table 3.


### 3.2.1. Initial Fracture Apertureb0

The pore-fluid pressure remained constant. The larger the initial aperture of the fracture was, the larger was the aperture and corresponding fracture permeability following deformation (Figure12).Figure 12
Fracture aperture and fracture permeability vs. initial fracture aperture.


### 3.2.2. Pore-Fluid Pressurep

As shown in Figure13, the initial fracture apertures were identical. Both, the fracture aperture and fracture permeability increased with increasing of fluid pressure.Figure 13
Fracture aperture and fracture permeability vs. pore-fluid pressure.


### 3.2.3. Initial Fracture Distanced0

The larger the fracture distance, the smaller the fracture density and the smaller the number of fractures in the same size unit. Therefore, when fluid pressure remained constant, the larger the fracture distance, the smaller the total apertures and the smaller the corresponding fracture permeability (Figure14).Figure 14
Fracture aperture and fracture permeability vs. initial fracture space.


### 3.2.4. Elastic Modulus

The elastic modulus of the matrix and fracture systems are related to each other. For the same seepage medium, it is necessary to discuss the change law for fracture permeability, when the two variables change simultaneously. As shown in Figure15, when the changes in effective stress are identical, the greater is the elastic modulus of the matrix system, and the smaller is the strain for the matrix system. Because the displacement values for the fracture system and the matrix system are equal in magnitude and opposite in direction, the deformation of the fracture system is smaller. Therefore, both fracture aperture and permeability decrease with increasing elastic modulus of the matrix system. When the elastic modulus of the fracture system increased, the change law did the same.Figure 15
Fracture aperture and fracture permeability vs. fracture modulus.


### 3.2.5. Effective-Stress Coefficient

There is also a specific relationship between the effective-stress coefficient of the matrix and the fracture system. However, there is no situation where one changes while the other remains unchanged. Therefore, for the same RE, it is necessary to discuss the change law of fracture permeability, when they change simultaneously. As shown in Figure16, the larger the effective-stress coefficient was, the smaller were the effective stress and the deformation of the matrix and fracture under the same total stress change and pore-fluid change. The larger the effective-stress coefficient was—following after the fracture has deformed—the smaller were the aperture and permeability.Figure 16
Fracture aperture and fracture permeability vs. effective-stress coefficient.


### 3.2.6. Fracture Angleβ

As the fracture angleβ increased, the total number of fractures in RE first increased and then decreased. Therefore, both the total fracture aperture and fracture permeability in isotropic media increase first and then decrease.Whenβ was 45°, the fracture permeability reached the maximum. When β was 0° and 90°, the fracture permeability was equal and reached the minimum (Figure 17).Figure 17
Fracture aperture and fracture permeability with fracture direction.


## 3.2.1. Initial Fracture Apertureb0

The pore-fluid pressure remained constant. The larger the initial aperture of the fracture was, the larger was the aperture and corresponding fracture permeability following deformation (Figure12).Figure 12
Fracture aperture and fracture permeability vs. initial fracture aperture.


## 3.2.2. Pore-Fluid Pressurep

As shown in Figure13, the initial fracture apertures were identical. Both, the fracture aperture and fracture permeability increased with increasing of fluid pressure.Figure 13
Fracture aperture and fracture permeability vs. pore-fluid pressure.


## 3.2.3. Initial Fracture Distanced0

The larger the fracture distance, the smaller the fracture density and the smaller the number of fractures in the same size unit. Therefore, when fluid pressure remained constant, the larger the fracture distance, the smaller the total apertures and the smaller the corresponding fracture permeability (Figure14).Figure 14
Fracture aperture and fracture permeability vs. initial fracture space.


## 3.2.4. Elastic Modulus

The elastic modulus of the matrix and fracture systems are related to each other. For the same seepage medium, it is necessary to discuss the change law for fracture permeability, when the two variables change simultaneously. As shown in Figure15, when the changes in effective stress are identical, the greater is the elastic modulus of the matrix system, and the smaller is the strain for the matrix system. Because the displacement values for the fracture system and the matrix system are equal in magnitude and opposite in direction, the deformation of the fracture system is smaller. Therefore, both fracture aperture and permeability decrease with increasing elastic modulus of the matrix system. When the elastic modulus of the fracture system increased, the change law did the same.Figure 15
Fracture aperture and fracture permeability vs. fracture modulus.


## 3.2.5. Effective-Stress Coefficient

There is also a specific relationship between the effective-stress coefficient of the matrix and the fracture system. However, there is no situation where one changes while the other remains unchanged. Therefore, for the same RE, it is necessary to discuss the change law of fracture permeability, when they change simultaneously. As shown in Figure16, the larger the effective-stress coefficient was, the smaller were the effective stress and the deformation of the matrix and fracture under the same total stress change and pore-fluid change. The larger the effective-stress coefficient was—following after the fracture has deformed—the smaller were the aperture and permeability.Figure 16
Fracture aperture and fracture permeability vs. effective-stress coefficient.


## 3.2.6. Fracture Angleβ

As the fracture angleβ increased, the total number of fractures in RE first increased and then decreased. Therefore, both the total fracture aperture and fracture permeability in isotropic media increase first and then decrease.Whenβ was 45°, the fracture permeability reached the maximum. When β was 0° and 90°, the fracture permeability was equal and reached the minimum (Figure 17).Figure 17
Fracture aperture and fracture permeability with fracture direction.


## 3.3. Influencing Factors for the Fracture Permeability Tensor

We introduced two rectangular coordinate systems (x′y′, xy) along the direction of fracture development and the direction of macroscopic pressure gradient, respectively. Assuming matrix permeability was km, the permeability tensor of RE in a fracture coordinate system is as follows:
(24)K′⟶=kf+km00km.After converting equation (24) to the macroscopic pressure-gradient coordinate system xy via coordinate transformation, the associated permeability tensor becomes [19, 38]
(25)K⟶=kxxkxykyxkyy=kf+2km2+kf2cos2βkf2sin2βkf2sin2βkf+2km2−kf2cos2β.Becausekm≪kf, equation (25) can be simplified to:
(26)K⟶=kfcos2βsinβcosβsinβcosβsin2β.The angle between the direction of the main permeability and thex-axis is as follows:
(27)β=12arctan2kxykxx−kxy+π2,β≤−π4,β=12arctan2kxykxx−kxy,−π4≤β≤π4,β=12arctan2kxykxx−kxy−π2,β≥π4.


### 3.3.1. Permeability-Tensor Componentskxx, kxy, and kyy Vary with Fracture Angle

The permeability-tensor componentkxx decreases with increasing fracture angle, kyy increases with increasing fracture angle, and kxy increases first but then decreases, with increasing fracture angle. When the fracture angle was 45°, the component kxy can reach its maximum value. When the fracture angle was 0° and 90°, the component kxy was 0 (Figure 18). It shows that as the fracture angle increased, the flow capacity in the x direction weakened gradually, while the fluid-flow capacity in the y direction gradually increased.Figure 18
Variation of the permeability tensor with fracture direction.


### 3.3.2. Variation Characteristic of Main Value Direction of Permeability

Within an isotropic seepage medium without fractures, the permeability was the same everywhere in the medium. According to equation (25), the permeability tensor ellipse of fractured media can be obtained (Figure 19). When a group of fractures develops in the media, the overall permeability of the media can show anisotropic characteristics. Moreover, the maximum main value direction for permeability was always along the fracture-development direction, while the minimum principal value direction was always perpendicular to the fracture-development direction.Figure 19
Anisotropic permeability tensor vs. fracture angle change.


### 3.3.3. Influence of Fracture Deformation on Permeability Tensor

Figure20(a) shows the variation of a permeability tensor vs. pore-fluid pressure in a seepage medium with a group of 0° fracture development. Figure 20(b) shows the variation of a permeability tensor with pore-fluid pressure in a seepage medium for a group of 45° fracture development. With increasing fluid pressure, the maximum principal value of medium permeability (length of elliptical semimajor axis) increases. However, the main value direction remains unchanged, and vice versa. This suggests that, given a seepage medium with a set of fracture development and isotropic matrix, the fracture deformation only alters the main value of the permeability.Figure 20
Anisotropic permeability tensor as a function of fluid pressure.
(a)
A set of 0° fractures(b)
A set of 45° fracturesUsing the two groups of fracture distributions shown in Figure21 as an example, when two groups of unsymmetrical fractures develop in RE (Figure 21(a)), the principal value of permeability increased with increasing fluid pressure. In addition, the direction of the permeability principal value gradually approached the group of fractures with strong pressure sensitivity (Figure 22(a)). Moreover, when two groups of fractures follow an orthogonal distribution (Figure 21(b)), or the sum of two groups of fracture angles was 90°, only the principal value changed (Figure 22(b)).Figure 21
Schematic diagram for two groups of fracture models.
(a)
0°~60°(b)
0°~90°Figure 22
Variation of main values of fracture permeability for two groups at different angles.
(a)
0°~60°(b)
0°~90°


## 3.3.1. Permeability-Tensor Componentskxx, kxy, and kyy Vary with Fracture Angle

The permeability-tensor componentkxx decreases with increasing fracture angle, kyy increases with increasing fracture angle, and kxy increases first but then decreases, with increasing fracture angle. When the fracture angle was 45°, the component kxy can reach its maximum value. When the fracture angle was 0° and 90°, the component kxy was 0 (Figure 18). It shows that as the fracture angle increased, the flow capacity in the x direction weakened gradually, while the fluid-flow capacity in the y direction gradually increased.Figure 18
Variation of the permeability tensor with fracture direction.


## 3.3.2. Variation Characteristic of Main Value Direction of Permeability

Within an isotropic seepage medium without fractures, the permeability was the same everywhere in the medium. According to equation (25), the permeability tensor ellipse of fractured media can be obtained (Figure 19). When a group of fractures develops in the media, the overall permeability of the media can show anisotropic characteristics. Moreover, the maximum main value direction for permeability was always along the fracture-development direction, while the minimum principal value direction was always perpendicular to the fracture-development direction.Figure 19
Anisotropic permeability tensor vs. fracture angle change.


## 3.3.3. Influence of Fracture Deformation on Permeability Tensor

Figure20(a) shows the variation of a permeability tensor vs. pore-fluid pressure in a seepage medium with a group of 0° fracture development. Figure 20(b) shows the variation of a permeability tensor with pore-fluid pressure in a seepage medium for a group of 45° fracture development. With increasing fluid pressure, the maximum principal value of medium permeability (length of elliptical semimajor axis) increases. However, the main value direction remains unchanged, and vice versa. This suggests that, given a seepage medium with a set of fracture development and isotropic matrix, the fracture deformation only alters the main value of the permeability.Figure 20
Anisotropic permeability tensor as a function of fluid pressure.
(a)
A set of 0° fractures(b)
A set of 45° fracturesUsing the two groups of fracture distributions shown in Figure21 as an example, when two groups of unsymmetrical fractures develop in RE (Figure 21(a)), the principal value of permeability increased with increasing fluid pressure. In addition, the direction of the permeability principal value gradually approached the group of fractures with strong pressure sensitivity (Figure 22(a)). Moreover, when two groups of fractures follow an orthogonal distribution (Figure 21(b)), or the sum of two groups of fracture angles was 90°, only the principal value changed (Figure 22(b)).Figure 21
Schematic diagram for two groups of fracture models.
(a)
0°~60°(b)
0°~90°Figure 22
Variation of main values of fracture permeability for two groups at different angles.
(a)
0°~60°(b)
0°~90°


## 3.4. Mechanism and Theoretical Explanation of Tensor Permeability Variation

According to the analysis results in Sections3.1, 3.2, and 3.3, under pressure sensitivity, not only will the principal value of the full tensor permeability of fractured anisotropic media change but the direction of principal permeability will also rotate. With increasing injection-production pressure difference or decreasing average pressure, the principal value direction will be inclined to the development direction of weakly sensitive fractures, and vice versa. The mechanism of this change is that several groups of fractures in the anisotropic strain medium produce nonlinear deformation on different degrees under pressure-sensitive conditions, which destroys the ratio relationship of the original permeability tensor.From the perspective of tensor theory, the permeability field is a symmetric tensor of second order. The change of a permeability tensor field caused by pressure sensitivity will inevitably lead to a “transformation effect” in the coordinate system of the permeability tensor. The transformation effect of the coordinate system has and only includes two aspects: (1) the expansion of the coordinate size, that is, the change of the principal value of permeability, and (2) the rotation of the coordinate system, that is, the change of the principal direction of permeability.


## 4. Conclusions

(1)
Under the constant-volume boundary, when fracture angle increased, the total fracture aperture and fracture permeability in the seepage medium first increased but subsequently decreased. Furthermore, the permeability of the 45° fracture system was the largest, and the others increased from a 0° to a 45° model, and then decreased from 90° to 45°(2)
Based on the analysis of the experimental results, in combination with the theory of seepage flow and rock mechanics, a novel analytical stress-dependent permeability model of fracture networks in different directions was derived. The calculation results of this model were in good agreement with the experimental results(3)
The novel analytical stress-dependent permeability model was used to study the influencing factors of the fracture aperture change. The study found the following: both total fracture aperture and the corresponding fracture permeability decrease with increasing fracture distance, elastic modulus of the matrix system, or effective-stress coefficient(4)
The influence of a single fracture and multiple fractures in different directions on the permeability tensor is analyzed by using the novel analytical stress-dependent permeability model. Our study indicates the following: for a seepage medium with only one set of fractures, the deformation of a fracture changing fluid pressure only affected the principal value of permeability. When two groups of unsymmetrical fractures developed in RE, the permeability principal value increased with increasing fluid pressure, and the permeability principal value direction gradually approached that of the fracture group with strong pressure sensitivity. When two groups of fractures were following an orthogonal distribution, only the principal value changed(5)
The pressure-sensitive equations of fractures in planes and longitudes throughout the study area are studied in this paper. However, there are both through fractures and nonthrough fractures in the actual reservoir. Therefore, the characterization method of pressure-sensitive permeability of nonthrough fractures can be further studied at a later stage



---

*Source: 7433547-2021-06-25.xml*