In fractured and cavity carbonate reservoirs, fractures form complex networks with multi-scale features, which include macro- and micro-fractures. It is notable that the fractures that significantly affect the fluid flow are those that connect caves because fault-related fractures have very high permeability potentials. Thus, fracture closure analysis focuses on fractures at this scale because of their important role in the conductivity of fractural and cavity carbonate reservoirs.

To clarify the calculation of the fracture closure pressure, the fracture is assumed to be a plate in the rock matrix, making the geo-stress field anisotropic. In detail, α is the dip angle of the fracture, and β is the angle between the strike of the fracture and the direction of maximum horizontal principal stress (Figure 1). Hence, the geo-stress field is a significant factor that affects the force analysis of the fracture plane. On the scale of the fracture surface, asperities are assumed to be hemispheroids with different radii and consistent spacing across the otherwise smooth and rigid surfaces. In particular, the asperities are independent, meaning that the deformation of a single asperity does not influence others even when the load increases until asperities on two opposite surfaces of the fracture come into contact.

With these assumptions, fracture closure occurs when opposing asperities come into contact under load, simultaneously causing variability in the fracture conductivity. As the load is applied, the distance between the asperities on the two plates decreases (Figure 2A), and the conductivity of the fracture depends on the fracture width only until the asperities on the two plates are in contact (Figure 2B). In these two stages, fractures still maintain considerable conductivity, so fracture closure slightly affects the production performance. After that, asperities become deformed, and the uncontacted fracture area is crucial for fracture conductivity (Figure 2C). In this situation, appropriate measures such as water injection are effective for preventing fractures from getting closer. Under a small load, asperities deform elastically, and the effective flow path decreases with the increasing load. When the load reaches a critical point at which the asperities yield, plastic deformation occurs along the fracture, and the fracture conductivity reaches the lowest value when the fracture closes completely (Figure 2D).

The onset of plastic deformation in materials is predicted by the yield criteria. Yield criteria are also called theories of yielding. Several yield criteria have been developed for ductile and brittle materials to predict the yielding behavior in simple ductile materials, and the Drucker–Prager failure criterion and Zienkiewicz–Pande failure criterion are fit for the material properties that may apply to fracture closure pressure calculation.

The theoretical differences between the Drucker–Prager failure criterion and the Zienkiewicz–Pande failure criterion lies in the mechanisms of yielding, which are caused by the hydrostatic pressure and strength differential effect. For the carbonate reservoir in the Shunbei Oilfield, the depth of the reservoir is approximately 8 km, and asperities on the fracture surface are significantly loaded by hydrostatic pressure, which means that the yielding caused by hydrostatic pressure is not negligible. In addition, the tensile strength of the rock is very low, approximately 0.1 times the compressive strength. Thus, rock materials are more likely to fail in tension than in compression (Aadnøy and Looyeh, 2019). This reveals that the strength differential effect is also important for fracture closure pressure estimation via the failure criterion. In contrast, the Zienkiewicz–Pande failure criterion is more suitable for fracture closure pressure estimation due to the yielding being caused by hydrostatic pressure and the strength differential effect.

The Zienkiewicz–Pande failure criterion was proposed based on a model for rocks and rock-like materials with multiple planes of weakness. The behavior of the assembly applies tensile and Mohr‒Coulomb shear limits on each such plane with possible strain dependence of the frictional properties. The model is applicable to the stability analysis of rock slopes. It is expressed as follows (Zienkiewicz and Pande, 1977; Li et al., 2017): f I 1 , J 2 = m 9 I 1 2 + n 3 I 1 − k + J 2 g 2 θ σ = 0 , (1) where g θ σ is the function of the stress Lode angle θ σ and is simply expressed as follows: g θ σ = 2 K ′ 1 + K ′ − 1 − K ′ sin ⁡ 3 θ σ . (2)

Parameters m , n , k , and K ′ are all functions of the internal friction angle ϕ 0 , which are expressed as follows: m = − tan 2 ⁡ ϕ 0 , (3) n = 2 c ¯ tan ⁡ ϕ 0 , (4) k = c ¯ 2 − a 2 ⁡ tan ⁡ ϕ 0 , (5) K ′ = 3 − sin ⁡ ϕ 0 3 + sin ⁡ ϕ 0 . (6) Here, a and c ¯ are defined as follows: c ¯ = 6 c ⁡ cos ⁡ ϕ 0 3 3 − sin ⁡ ϕ 0 , (7) a = c ¯ 24 ⁡ tan ⁡ ϕ 0 . (8)

As mentioned above, since fracture aperture is sensitive to effective load, force analysis of fluids in fractures is fundamental to studying fracture closure via stress variation. In our research, the geo-stress field and fluid pressure are introduced in the force analysis.

Derived from the fracture model we built in Figure 1, the effective normal stress on fracture σ n e is as follows: σ n e = σ v ⁡ cos ⁡ α + sin ⁡ α σ H ⁡ sin ⁡ β + σ h ⁡ cos ⁡ β − P f , (9) where σ v is the gravity of the overlying rock, which is determined by the density of the overlying rock and its depth. σ h and σ H are the minimum horizontal principal stress and the maximum horizontal principal stress, respectively, and these two stresses depend on the tectonic stress and horizontal stress derived from the vertical principal stress. P f is the fluid pressure, which is the force on the fluid in the fracture.

The limit of yielding is assumed to be equal to the effective normal stress and vertical principal stress: σ s = σ n e = σ z . (10)

In subsurface formations that are not subjected to significant tectonic forces, σ 1 is in the vertical direction due to the lithostatic pressure of the overburden, but significant compressional forces can result in σ 1 being oriented in or close to the horizontal plane. In our case, σ v = σ 1 , σ H = σ 2 , and σ h = σ 3 (Rackley, 2017).

Additionally, triaxial stresses are expressed as follows: σ x = σ y = μ 1 − μ σ z , (11) where μ is the Poisson’s ratio.

Then, the Zienkiewicz–Pande failure criterion is used in closure pressure analysis. From Equations 1–11, the closure pressure of the fracture derived from the Zienkiewicz–Pande failure criterion is expressed as follows: P f = σ v ⁡ cos ⁡ α + sin ⁡ α σ H ⁡ sin ⁡ β + σ h ⁡ cos ⁡ β − N 2 − 4 M K − N 2 M , (12) Parameters M, N and K in Equation 11 are expressed as Equations 12–15. M = 4 3 ⁢ sin 2 ⁡ ϕ 0 3 − sin ⁡ ϕ 0 2 1 + μ 1 − μ ⁢ 2 + 2 μ − 1 1 − μ ⁢ 36 − 34 ⁡ sin ⁡ ϕ 0 ⁡ sin 2 ⁡ 3 θ σ − 4 ⁡ sin 2 ⁡ ϕ 0 ⁡ sin 2 ⁡ 3 θ σ 96 + 12 ⁡ cos 2 ⁡ ϕ 0 , (13) N = 8 c ⁡ sin ⁡ ϕ 0 ⁡ cos ⁡ ϕ 0 3 − sin ⁡ ϕ 0 2 , (14) and K = 12 c 2 ⁡ cos 2 ⁡ ϕ 0 1 − 576 ⁡ tan ⁡ ϕ 0 576 ⁡ tan ⁡ ϕ 0 3 − sin ⁡ ϕ 0 2 . (15)

Based on the fracture closure pressure equation, the pressure is determined by three groups of parameters: the geo-stress field parameters, which include the gravity of the overlying rock, the minimum horizontal principal stress, and the maximum horizontal principal stress; the rock mechanics parameters, which involve Poisson’s ratio, the internal friction angle, and rock cohesion; and the fracture parameters, such as the dip angle of fracture and the angle between the fracture strike and the maximum principal stress. Higher geo-stress values result in higher fracture closure pressure; thus, we focus on the influence analysis of rock mechanics and fracture parameters.

For the oil production wells in the fractured-cavity reservoir, fractures are the dominant factor in production performance, and the aperture of the fracture determines its conductivity from the reservoir to the wellbore. In our research, oil production wells located in the same fault and exploitation unit were selected as typical wells to evaluate the fracture closure pressure (Figure 5).

From geophysical data, σ H , σ h , and σ v are 180 MPa, 130 MPa, and 150 MPa, respectively, for the Ordovician carbonate reservoir in the Shunbei Oilfield, and the direction of maximum principal stress is N54°E about the fault zone. Rock mechanics parameters were obtained by experiments and sonic log interpretation. The dip and strike of fractures were derived from well log and seismic interpretation images. The values of fracture parameters for fracture closure pressure are listed in Table 3.

Finally, we obtained the fracture closure pressures for these three wells, which are 69.28 MPa, 61.91 MPa, and 73.11 MPa for wells S1, S2, and S3, respectively (Figure 6). It is observed that the fracture closure pressure of S3 is the highest. The reason, inferred from the analysis of the influencing factors, is that the fractures of well S3 have the largest angle between the fracture strike and the maximum principal stress. From the perspective of tectonic genesis, this may be caused by the stress transition zone. Correspondingly, the fracture closure pressures of wells S1 and S2 are lower due to smaller angles. The evaluation of the influencing factors reveals that a large dip angle leads to a lower closure pressure. However, for well S3, although the dip angle of the fracture is the smallest, the fracture closure pressure is not as high as that of other wells. This is because the strike of the fractures is almost parallel to the direction of the maximum principal stress. Based on the pressure calculation, the difference in fracture closure pressure among the wells located in this fault in the Shunbei Oilfield is mainly determined by the strike of the fractures.

Since we estimated the closure pressure for the Shunbei Oilfield based on the Zienkiewicz–Pande failure criterion, stress field parameters, and rock mechanics parameters of the Shunbei Oilfield, the production performance of the oil wells after fracture closure was also discussed to verify the closure pressure estimate.