Well productivity evaluation in deformable single-fracture media

Evaluations and interpretations of reservoir productivity are frequent in geothermal, groundwater and hydro- carbon research and applications. In this study, we consider the closure of fractures around production wells due to compaction that can affect the productivity value, i.e. the ability of subsurface formations for transporting the desired fluid to a borehole. We introduce analytical tools to evaluate and predict changes in the productivity of a deformable fractured porous media. We propose analytical models for three geometries: a rectangular fracture with zero/non-zero orientation and a circular fracture with zero orientation to maximum horizontal stress. An advanced numerical model is utilised to evaluate the impact of spatial variation of fracture aperture induced by the fracture deformation on well productivity. The developed analytical solutions using a uniform fracture aperture always either over- or un- derestimate the production rate. Hence, an equivalent aperture model is developed for the fracture with aperture distribution under variable contact stresses to circumvent this problem. The proposed equivalent aperture model reduces the average and maximum errors of production-rate pre-diction from 28% to 0.6% and from 116% to 25%, respectively. We further employ the proposed model for sensitivity analyses to illustrate the impacts of in-situ and human-controlled parameters on productivity reduction. These analyses present that the interactions among initial reservoir pressure, fracture orientation, fracture stiffness, and well pressure control productivity reduction behaviours and the maximum productivity reduction values.


Introduction
Many production wells in geothermal, groundwater, and hydrocarbon reservoirs suffer from productivity reduction during depletion (Ren and Guo, 2018;Kadeethum et al., 2019b). The combined hydromechanical effect may result in this reduction because when fluid is withdrawn from the system, the effective stress increases, which in turn, reduces the fractures or rock matrix conductivity (Kadeethum et al., 2018b(Kadeethum et al., , 2019a. In the case of highly fractured reservoirs, the closure of natural fractures around the production wells due to the hydromechanical effect is expected to be one of the main reasons for the productivity reduction (Kadeethum et al., 2019b). This productivity reduction may significantly reduce the energy/water production capability (Wang et al., 2017). In some cases, stimulation treatments such as hydraulic fracturing or acid injection are required to maintain productivity and sustain heat or hydrocarbon production (Legarth et al., 2003).
Different numerical modelling approaches have been developed to study coupled flow and solid deformation in fractured porous media (Matthäi et al., 2010;Castelletto et al., 2016;Vasilyeva et al., 2019a,b;Vik et al., 2018;Adler et al., 2013). The interplay between flow and rock deformation can significantly influence the flow and transport properties of subsurface systems Vik et al., 2018;Bisdom et al., 2016). For example, Salimzadeh et al. (2018c) and  illustrated neglecting this coupling results in large errors for predicting the lifetime and the net energy production of fractured geothermal systems. While these methods provide flexibility in capturing detailed geological features, they suffer from high computational time and cost. Analytical solutions become an alternative method as they demand less computational time and resources (Streltsova, 1987;Sedghi et al., 2018;Dewandel et al., 2018;Adler et al., 2013). Past endeavours of the analytical solutions used to predict well productivity in fractured porous media involve both the continuum (Horne, 1995;Bogdanov et al., 2003) and discrete fracture-matrix approaches (Ibrahim et al., 2006;Bello et al., 2010;Kanfar et al., 2017). So far, however, there are not many published analytical solutions that https://doi.org/10.1016/j.geothermics.2020.101839 Received 31 July 2019; Received in revised form 11 March 2020; Accepted 16 March 2020 represent the coupled hydromechanical effect on well productivity in fractured porous media considering both fractures and the rock matrix explicitly. To identify the parameters controlling the productivity, and subsequently provide mitigation solutions fast mathematical models combining hydraulic and mechanical processes are needed. The previous analytical models that include the effect of stress-induced fracture conductivity only provide solutions in a rectangular domain with a fracture orthogonal to minimum horizontal stress (Tabatabaie et al., 2015). These models, however, do not account for the interaction between fluid pressure and fracture aperture (Tabatabaie et al., 2017). This interaction leads to a spatial variation of fracture aperture, which can significantly influence fluid flow behaviour and well productivity (Kadeethum et al., 2018a;Salimzadeh and Khalili, 2015).
To this end, we aim to develop and extend both analytical and numerical solutions, coupling the non-linear relationship between fracture pressure and aperture, and describe the advantages/disadvantages of each model. The proposed models can serve as a tool to evaluate productivity reduction in fractured porous media. This paper is structured as follows. First, the derivations of three analytical solutions and numerical model. Next a verification of the developed analytical solutions against the robust numerical model is presented. Then, an equivalent aperture model that can capture the preferable trait of the numerical method and provide a competitive component for the analytical solutions is introduced and tested. Furthermore, the factors controlling the productivity reduction are identified.

Methodology
In this section, we develop analytical solutions for three geometries (see Fig. 1) for single-phase flow at steady-state condition assuming the fluid flow is in a bi-linear flow regime. This implies that the fluid flows from domain boundaries to the fracture through the rock matrix in ydirection; subsequently, the fluid flows along the fracture to the well in x-direction. We also introduce the numerical model used to study the impact of solid deformation on fracture aperture and fracture flow.

Analytical solutions
This section elaborates governing equations and their solutions, which are utilised to predict production rate and pressure at each particular point in the fracture. Firstly, solutions for a rectangular fracture in a rectangular domain for each specified orientation are presented. Subsequently, a solution of a circular fracture located in the cylinder domain with 0°orientation is derived.

Rectangular fracture with zero orientation
The geometry of a rectangular fracture with 0°orientation (respect to x-direction) is shown in Fig. 1a and b. According to a model symmetry, only 1/8 of the model is required to be solved (the section is shown by red colour in Fig. 1a and b). As the model is at a steady-state regime, Darcy's law in the fracture space is utilised using the following equation: where q f (x) is fracture flow rate at any given x, μ is fluid viscosity, p f (x) is the fracture pressure at any given x, k f and A are fracture permeability and cross-sectional area, respectively, and read as follows: where a cal is the calculated fracture aperture, and h is the fracture height. In this study, fractures are in contact, which means two sides of fracture are physically touching each other. Thus, fracture asperities are the primary media to provide flow channels that can be represented as an average contact aperture (Jaeger et al., 2009;Bisdom et al., 2016).
To determine this value, the empirical non-linear relationship, known as the Barton-Bandis model, is used: (Bandis et al., 1983;Barton et al., 1985).
where a 0 is the fracture aperture at zero contact stress, which can be approximated by = a a b 0 assuming there is no residual aperture when the contact stress goes to infinity. The a and b are model parameters, and σ′ n is a normal component of the contact stress, which is defined as: where σ n is a normal component of the far-field stresses that act on a fracture surface. From Fig. 1, all of the matrix flow is collected through a single fracture connected to the production well. The fluid flow in fracture at any given x can also be written as follows: where x f is the fracture half-length, q m (x) is matrix flow rate at any given x. Subsequently, flow in the matrix is also governed by Darcy's law as: where k m is the matrix permeability, p i is the matrix pressure at the boundary (y = l m ), l m is the distance from the boundary to the fracture. The following variable is used: and one constant is formed as: The boundary condition at x = 0 represents the differential pressure between the pressure at the boundary and the well pressure (p f (x = 0)): and the p at x = x f reflects a no-flow boundary condition: By enforcing the boundary conditions (Eqs. (11) and (12)) in Eq. (10), we then seek a solution of p(x) as: e e x x f x f x f x x f Combining Eqs. (1) and (13), then applying it for the whole domain, the well production rate is read: where q total is the total production rate.

Rectangular fracture with non-zero orientation
In the case that fracture has the angle of θ with respect to x-direction (0 < θ < 90) as presented in Fig. 1c and d, governing equations are adapted from the previous section to be in compliance with model geometry, i.e. the fracture is located in x′-coordinate instead of x-coordinate. Eqs. (1), (6), (8, 11), and (12) are adjusted to x′-coordinate, and only 1/4 of the system (shown in Fig. 1c and d in red) is solved as follows: The normal component of the far-field stresses with respect to the fracture plane, σ n , in Eq. (5) is determined through the following equation: where σ 1 and σ 2 are the minimum and maximum horizontal stresses (yand x-direction), respectively. The constant ψ is also modified as: The governing equation for a fracture with θ°orientation is summarised as follows: Again, by enforcing Eqs. (11) and (12) to Eq. (18); and subsequently, utilising Lie symmetry analysis (Hereman, 1977;Liu et al., 2009), we then seek a solution of p x ( ) as follows: )) f where is the Gaussian or ordinary hypergeometric function, and the definition of coefficients c 1 to c 4 and z 1 to z 12 can be found in Appendix A. Using Eqs. (1) and (19) the production rate for the domain shown in Fig. 1c and d is read:

Circular fracture with zero orientation
The analytical solution for the flow into a circular fracture in a cylindrical domain is given. The origin of the coordinate is in the middle of the circular fracture. According to the symmetry of the model, only 1/8 of the domain is used (see Fig. 1). Flow is steady-state; hence, Darcy's law in the fracture is written as: where q f (r) is fracture flow rate at any given r and p f (r) is fracture pressure at any given r. Since in this case, fracture geometry is circular, the fracture cross-sectional area is calculated as follows: From Fig. 1e and f, all of the matrix flow is accumulated by a single circular fracture. Therefore, fluid flow in fracture at any given r can also be written as shown below: where R is fracture radius, q m (r) is matrix flow rate at any given r in a rdirection. The flow in the matrix is steady-state expressed by: We modify Eqs. (8) and (9) to solve this system of equations, and the final ordinary differential equation is presented as follows: Two boundary conditions, the first one is a constant well pressure as presented below: where r 1 is the well radius and the second one is the no-flow boundary as: By enforcing Eqs. (26) and (27) to Eq. (25); and subsequently, utilising Lie symmetry analysis (Hereman, 1977;Liu et al., 2009) to solve this ordinary differential equation results in: where X n (o) and Y n (o), n = 0, 1, are Bessel function of the first and the second kind, respectively, and coefficients o 1 to o 3 are described in Appendix A (Abramowitz and Stegun, 1964;Finch, 2003). The production rate for the whole domain reads: T. Kadeethum, et al. Geothermics 87 (2020) 101839 Note that this set of equations is valid when 0 < r 1 < R.

Numerical model
The analytical solutions discussed in the previous section do not have the capability to consider the fracture aperture spatial variation.
We utilise a hydromechanical (HM) numerical model developed by Salimzadeh et al. (Salimzadeh et al., 2017a, 2017b to evaluate well productivity reduction resulting from fracture aperture closure. In this model, fractures and wells are modelled as high-permeability surfaces (2D) and lines (1D) in the 3-Dimensional matrix domain, respectively.
Fractures can be modelled as either internal walls or split surfaces; however, when fracture surfaces are in contact and propagation of fractures is not likely occurred, the internal wall method is preferred because it is computationally cheaper. The flow field is continuous across fracture and rock matrix, and it is written as: where k m is the matrix permeability tensor, ρ is the fluid density, g is a vector of gravitational acceleration, α is the Biot coefficient, K is rock bulk modulus, K s is solid bulk modulus, φ is rock porosity, c f is fluid  .where k w is well permeability and r 1 is well radius. To mimic mechanical deformation due to the changing pressure and represent aperture variation due to effective normal stress, the following equation is used: where D is a drained stiffness matrix, ε is the strain, α is the Biot's coefficient, σ′ is the stress tensor, n is a unit vector perpendicular to the external boundaries and I is the second-order identity tensor. Eqs. (30) and (31) are discretised based on continuos Galerkin finite element method for spatial domain. The temporal space is discretised by the backward Euler scheme. The primary variables are displacement vector, u, matrix pressure, p m , and fracture pressure, p f . The model is implemented in Complex Systems Modelling Platform (CSMP), an object-oriented application programme interface (Matthäi et al., 2010;Nick and Matthäi, 2011). Algebraic Multigrid Methods for Systems (SAMG) is utilised to solve the system of equations in each time-step (Stüben et al., 2017). Quadratic elements, second-order approximation, are used in all cases throughout this study.

Fracture aperture representation
Two fracture aperture representations are used in this study: uniform aperture and deformable (non-uniform) fracture aperture. The uniform aperture model uses the fracture aperture calculated from constant pressure value, initial reservoir pressure (p i ) or well pressure (p wf ), which are named a pi or a pwf , respectively. For the deformable case, the fracture aperture at each point is dynamically computed from the corresponding contact stress value, which is a function of p f and u at that point . The deformable aperture varies spatially, and examples of deformable fracture are illustrated in Fig. 2 for both rectangular and circular fracture geometries (p wf = 500 psi). The rectangular fracture has an aperture variation in two directions (horizontal and vertical directions), while the circular fracture has the aperture variation only in one direction (the radial direction). This observation can be explained as the well geometry of the circular fracture intersects the fracture as a point at the centre. In contrast, the well in the rectangular fracture intersects the fracture as a line through Fig. 4. Deformable aperture effect on reservoir deliverability: rectangular fracture with (a) 0°, (b) 30°, (c) 45°, and (d) 60°orientation (θ) model, the deformable aperture results remain between two uniform aperture results, which calculate fracture aperture as a function of the initial reservoir pressure and the well pressure.
T. Kadeethum, et al. Geothermics 87 (2020) 101839 a centre axis. The summary of the discussed representations is shown in Table 1. In short, a pi and a pwf models lead to a cal and k f becoming a function of p i or p wf . The deformable aperture, on the other hand, captures the fractures that response to the coupled hydromechanical effect. Therefore, the a cal and k f become the function of p f and u at each specific point. Note that a pi and a pwf are also referred to the a cal that is calculated using p i or p wf value, respectively. For the analytical solutions, only the uniform aperture calculation method is employed, while both representations are used in the numerical simulations.

Productivity index evaluation
The productivity index (J), which shows the ability of subsurface formations to deliver a desirable fluid through a borehole, is defined as (Porges, 2006): where p e is external boundary radius pressure and p wf is well pressure. The production rate and well pressure are straightforward to measure and determine; on the contrary, external boundary radius pressure requires a more complex method to estimate. Therefore, our problem is simplified by representing the system for a steady-state condition and calculating J using constant pressure at the boundaries, p i . Subsequently, Eq. (32) is adjusted as follows:

Table 3
Fracture aperture comparison: p wf = 500 and 2000 psi for deformable cases; aperture (p i ) represents fracture aperture that is calculated using the reservoir initial pressure while fracture aperture in deformable aperture model is computed using fracture pressure. T. Kadeethum, et al. Geothermics 87 (2020) 101839 where Δp is the difference between p i and p wf or in this case p f (0). We calculate the J of the production well by following procedures; (i) we prescribe the pressure value, p i , at the boundaries shown in Fig. 1 while we assume no-flow boundary conditions at the rest of the domain boundaries; (ii) well pressure, p wf , is strongly enforced at the well (see green circle in Fig. 1); (iii) the production rate corresponding to the specified p wf is determined, and (iv) J is calculated through Eq. (33). By keeping p wf constant and changing boundaries pressure (p i ), the model represents a pressure depletion through fluid production. On the other hand, by retaining p i and altering p wf , the evolution of J with different drawdowns (p i − p wf ) can be observed. This process ensures that J is calculated consistently. Moreover, the coupled hydromechanical effect on reservoirs deliverability (fracture aperture alteration) can be identified.

Results and discussion
The outline of this section is summarised as follows: firstly, the verification of the analytical solutions with the established numerical model is presented. Secondly, the results of uniform aperture and deformable aperture cases are compared and discussed. Thirdly, the equivalent aperture model is developed to imitate the deformable aperture effect. Fourthly, the equivalent aperture model is tested with different variable ranges. The first range uses input parameters that are inside a development dataset -hereafter referred to as interpolation. While the second one uses input parameters that are outside the development dataset -hereafter referred to as extrapolation (see Table 4). Finally, a sensitivity analysis using the equivalent aperture model is performed to evaluate the impact of initial reservoir pressure, fracture orientation, fracture stiffness model parameters, and well pressure on system productivity reduction.

Analytical solution verification
We consider two cases: the first case calculates aperture as a function of the initial reservoir pressure, a pi . At the same time, the other determines fracture aperture as a function of the well pressure, a pwf as discussed in the previous section and referred to as the uniform aperture representation. Input parameters for analytical and numerical models are presented in Table 2. In this verification, four p wf values are applied to the rectangular fracture model, and five p wf values are applied to the circular fracture model.
Eqs. (14) and (20) are used to calculate production rate for the rectangular fracture model with 0°and θ°orientation, respectively. Three fracture orientations, 0°, 30°, and 60°, are considered in this verification. Fig. 3a-c show good agreements between the numerical and analytical production rate results. The same agreements are illustrated in Fig. 3d for the circular fracture model using Eq. (29).  6. (a) Deformable aperture impact on reservoir deliverability: circular fracture with 0°orientation model, the deformable aperture results remain between two analytical solutions, which calculate fracture aperture as a function of the initial reservoir pressure and the well pressure. (b) fracture aperture variation with different well pressures. Aperture at initial reservoir pressure is 0.233 mm, and aperture at the well pressure for each well pressure case is located at lower bound of that case.

Deformable vs. uniform aperture representations
Input data of this analysis is presented in Table 2. Rectangular fracture geometry is simulated for four orientations, 0°, 30°, 45°, and 60°, and each geometry has four different well pressures, p wf , presented in Fig. 4. There are two key results from this investigation: (1) the relationship between p wf and production rate in the deformable aperture model is non-linear; (2) all deformable model production rate results are located between a pi and a pwf models with a tendency towards the lower bound (a pwf model).
This behaviour can be explained by the aperture distributions, as shown in Fig. 5. This figure uses Kernel density estimation normalised by the maximum frequency in each case. As observed in Fig. 5, the fracture apertures of a pi model are much greater than the maximum aperture of the deformable model. Hence, the production results between deformable and a pi models are more diverse than the results between deformable and a pwf models as demonstrated in Fig. 5 and Table 3. The results also show that when p wf is reduced, the variation of the aperture is increased. Note that apertures of a pwf model are the lower bounds of deformable apertures for each particular p wf case because the lowest pressure point in the system is at the well.
The production rate and fracture aperture results of the circular fracture geometry for seven p wf values are presented in Fig. 6. Two observations noted for the rectangular fracture geometry are held. The further observations can be made: (1) fracture apertures change dramatically near the well and evolve gradually far from the well. This means that the near-well domain may dictate the system behaviour and need more investigation to understand this behaviour. (2) The differences between the maximum and minimum apertures are increased when p wf is increased. Hence, reducing p wf increases the differences between the deformable fracture and a pwf model production rate results.
Sets of numerical experiments have been performed with input parameter ranges shown in Table 4 -interpolation to elaborate on the effect of the deformable aperture model. The differences between production rate results of deformable aperture model and uniform aperture model, calculated at p i and p wf , are shown in Fig. 7. The production rates calculated by a pi model are always greater than those of the deformable model. The production rates calculated by a pwf model, on the other hand, are always less than those of the deformable model. The   Fig. 7. (a) the production rate (q) of a pi and a pwf models vs. q of deformable fracture model and (b) q difference distribution; average errors of a pi and a pwf models are 31% and 10%, respectively; maximum errors of a pi and a pwf models are 113% and 53%, respectively; (c) a pi and a pwf values vs. equivalent apertures and (d) aperture difference, a contrast between equivalent apertures and a pi or a pwf values, distribution; average errors of a pi and a pwf models are 29% and 8%, respectively; maximum errors of a pi and a pwf models are 154% and 24%, respectively.
results of a pwf model are closer to the results of the deformable model than those of a pi model. Errors of a pi and a pwf models are calculated as follows: where q is the production rate. a pwf and a pi models have the average error of 10% and 31%, respectively, and the maximum error of 53% and 113%, respectively. a pwf model is more accurate than a pi model because a pwf model takes the effect of the fracture aperture reduction due to changing of fracture pressure into account.

Equivalent aperture model development
We propose an equivalent aperture model that captures the behaviour of the deformable aperture model. The procedure starts with randomly generating sets of input parameters within the ranges shown in Table 4 -interpolation. Then, the deformable aperture model is simulated, and its production rate is obtained. The equivalent aperture that produces the same production rate as the deformable aperture model is achieved through an iterative process using the analytical solutions. This procedure begins with an initial guess of fracture aperture value and keeps changing the aperture value until the tolerance criterion is met. We define the tolerance criterion as a difference between the production rate of the deformable aperture model and the calculated production rate using the analytical models. The equivalent aperture model is established by multivariable regression and Box-Cox transformation.
The multivariable linear regression of a eq contains four independent Fig. 8. Model verification -input data is inside the development range: (a) equivalent aperture model (b) a pi and a pwf models; average errors of equivalent aperture, a pi , and a pwf models are 0.7%, 28%, and 7%, respectively, and maximum errors are 15%, 80%, and 24%, respectively; input data is outside the development range: (c) equivalent aperture model (d) a pi and a pwf models; average errors of equivalent aperture, a pi , and a pwf models are 0.6%, 25%, and 6%, respectively, and maximum errors are 25%, 116%, and 46%, respectively. q is production rate.  Kadeethum, et al. Geothermics 87 (2020) 101839 Fig. 9. Productivity index, J, reduction sensitivity analysis; this figure illustrates the calculation procedure of productivity, J, reduction: (a) well pressure vs. production rate plot -grey arrows represent the production rate different between a pi and deformable aperture models, (b) well pressure vs. J reduction. Fig. 10. J reduction of (a) rectangular fracture 0°geometry, (b) rectangular fracture 60°geometry, and (c) circular fracture 0°geometry.
T. Kadeethum, et al. Geothermics 87 (2020) 101839 variables, and the equation is read: where c 1 = 3.645625, c 2 = 0.115169, c 3 = 0.934210, c 4 = 2.895228, c 5 =−0.141167 and c 6 = 0.6616427. Note that we have selected this dimensionless form as the most accurate one among multiple tested relationships. We use R, an open-source software environment for statistical computing and graphics, to build this regression model (Kleinbaum et al., 2007;R Core Team, 2017;Levine et al., 2001). More details of a eq development, R-squared, and residual analysis are described in Supplementary Information.

Equivalent aperture model verification
To verify the proposed model (Eq. (35)) against scenarios that are not considered in the model development an extra set of test dataset (Table 4 -extrapolation) are utilised. Four rectangular fracture geometries, 0°, 30°, 45°, and 60°orientation, and one circular fracture geometry, 0°orientation, are used in this verification.
The results of 356 simulations of the interpolation range (Table 4) are presented in Fig. 8a-b, and the results of 704 simulations of the extrapolation range (Table 4) are presented in Fig. 8c-d. Moreover, the summary of the R-squared, average, and maximum errors of a eq , a pi , and a pwf models is presented in Table 5. These results demonstrate the a eq model superiority over the a pi and a pwf models. Furthermore, the a eq model provides an acceptable predictive performance.

Sensitivity analysis of productivity index reduction
After the a eq model is established this model is utilised to present impacts of p i , θ, and a 0 , on system productivity (J) behaviour through a sensitivity analysis. Each parameter is investigated by interacting with different p wf to evaluate the impact of the interaction between in-situ and human-controlled variables. These variables are selected because the p i directly influences the fracture aperture value. Moreover, θ influences the normal contact stress, which in turn affects the fracture aperture. Lastly, a 0 not only defines the initial fracture aperture but also dictates the deformable behaviour.
The sensitivity analysis is performed based on one-factor-at-a-time method (Daniel, 1973). Input parameters for the base case are presented in Table 2. For each case, p wf is varied from 0.1 to 0.9 of its p i value. Eqs. (14), (20), (29) and (35) are used to calculate the production rate depending on the model geometry. After production rate is obtained from the mentioned equations, J is calculated using Eq. (32) for two cases. The first case production rate is obtained from the a pi model, and the second one uses the a eq model. Subsequently, J reduction is read as follows: Eq. (36) presents the productivity reduction as a result of the effective stress dependency of the fracture aperture; in other words, the coupled hydromechanical effect on the fracture conductivity. To elaborate this effect, Fig. 9a illustrates a linear relationship between p wf and production rate of a pi model but a non-linear relationship between p wf and production rate of a eq model. Thus, the coupled hydromechanical effect makes the relationship between well production rate, q, and p wf nonlinear. Moreover, the J reduction has non-linear relationship with p wf as presented in Fig. 9b. The effect of interaction between p i and p wf on J reduction is presented in Fig. 10. Note that p i is varied from 0.1 to 0.9 of the minimum horizontal stress, σ 1 . Fig. 10 illustrates that the maximum J reduction in all models occurs when p wf is the lowest, but p i is the highest. The magnitude of J reduction between rectangular and circular fracture models are very different. Since the circular fracture model has a single point of production in the middle, it suffers from the fracture aperture reduction that acts as a bottleneck and prevents the flow from the fracture boundary to the production point. Although the rectangular fracture model is affected by the same phenomenon, it contains a line producer, which provides more slots to produce fluid from the reservoir.
The effect of interaction between θ and p wf on J reduction is presented in Fig. 11a for the base case; θ is varied in the range of 5°to 85°. The high J reduction zone is located where θ is approximately 45°(σ 1 and σ 2 influence σ n equally.), and p wf is the lowest. This is because the difference between the initial aperture and the deformed aperture becomes the largest for the lowest value of p wf . Moreover, when θ is increased, σ 2 impact becomes larger and causes the initial J smaller. As a result, the difference between the initial J and the deformed J is small.
Next we investigate the effect of interaction between θ and far-field stresses, σ 1 and σ 2 , on J reduction. σ 1 and σ 2 are 1.4 of their base case values ( Table 2). The results are presented in Fig. 11b. The high J reduction zone is located at the same location at its base case. Moreover, T. Kadeethum, et al. the contour lines are also similar. However, the highest J reduction value of this case is higher than that of the base case because of increasing normal contact stress due to the higher σ 1 and σ 2 values.
The effect of interaction between a 0 and p wf on J reduction is presented in Fig. 12; a 0 is varied in the range of 0.5 to 1.5 of its base case value (Table 2). Since a 0 is estimated using a 0 = a/b there is no unique pair of a and b that gives a particular a 0 ; therefore, the investigation is divided into two parts. The first part, a value is varied while b is fixed. The second part, on the other hand, a is fixed, but b is varied. Note that both cases have the same a 0 . There are two main observations from this investigation: (1) the maximum J reduction occurs when a 0 and p wf are smallest, because a p is the smallest when a 0 and p wf are at the minimum Fig. 12. J reduction of the rectangular 0°orientation fracture model; (a) a is varied, but b is fixed and (b) a is fixed, but b is varied; the rectangular 60°orientation fracture model; (c) a is varied, but b is fixed and (d) a is fixed, but b is varied; the circular 0°orientation fracture model; (e) a is varied, but b is fixed and (f) a is fixed, but b is varied.