# Pressure transient analysis for stress-sensitive fractured wells with fracture face damage

Ruben Figueroa Hernandez, Anis Younes, Marwan Fahs, Hussein Hoteit

► To cite this version:

Ruben Figueroa Hernandez, Anis Younes, Marwan Fahs, Hussein Hoteit. Pressure transient analysis for stress-sensitive fractured wells with fracture face damage. Geoenergy Science and Engineering, 2023, 221, pp.211406. 10.1016/j.geoen.2022.211406. hal-04296285

HAL Id: hal-04296285

https://hal.science/hal-04296285v1

Submitted on 24 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Pressure transient analysis for stress-sensitive fractured wells with fracture face damage

Ruben Figueroa Hernandez<sup>a</sup>, Anis Younes<sup>b</sup>, Marwan Fahs<sup>b</sup>, Hussein Hoteit<sup>a, *</sup>

## ARTICLE INFO

**Keywords:**

Pressure transient analysis

Well testing

Hydraulic fractures

Fracture closure

Fracture face skin

Formation damage

Stress-dependent permeability

Type-curve

## ABSTRACT

Fracture networks improve wells' deliverability of fluids from underground reservoirs. However, conductive fractures may lose efficiency because of fracture closure or formation damage that needs to be monitored and detected for proper reservoir management and intervention. Fracture hydraulic conductivity is a function of fracture geometry and its hydraulic aperture, which tends to degrade during operations because of fracture compaction and closure. Furthermore, pore-clogging of fracture walls (faces) results in a damaging skin zone that may reduce fracture-matrix transmissibility. Current Pressure Transient Analysis (PTA) models in the literature do not consider the interconnected combined effect of stress-dependent conductivity of fractures and fracture face skin. Modeling fracture closure without accounting for the skin effect may mislead the PTA interpretations of well performance. This work proposes a rigorous PTA model to account for fracture closure combined with fracture face skin for single and multistage hydraulic fractures with uniform and complex geometries. The solution approach is based on subdividing the domain into three flow regions, matrix, damage zone, and hydraulic fracture (HF). The proposed model is verified extensively using a commercial simulator and other less comprehensive models from the literature. The pressure response in fractured wells with stress-sensitive fractures is analyzed at early, middle, and late times. In each flow regime, we identify pressure signals to detect fracture closure by incorporating fracture closure and skin effects. Results show that the fracture face skin has a significant effect on fracture performance as it introduces an additional pressure drop that triggers an earlier and more pronounced occurrence of fracture closure. We proposed a semi-log approach to identify fracture closure for low compressible fractures and a pseudo-radial simplification to generate late-time response curves. The proposed method serves as a PTA diagnostic tool to detect fracture conductivity degradation due to fracture closure and formation damage.

## 1. Introduction

Hydrocarbon production from tight reservoirs requires specialized techniques, including horizontal drilling and Multistage Hydraulic Fractured Wells (MHFW). Hydraulic fracturing consists of a series of multistage high permeability fractures created by the injection of high-pressure fluids, followed by particle (proppant) placement to prevent fracture closure due to stresses exerted by the surrounding formation (Cardona et al., 2021; Fan et al., 2011; King, 2012; Male et al., 2014; Patzek et al., 2013). In other situations, hydraulic fractures could be induced unintentionally due to uncontrolled fluid injection in secondary recovery methods such as water flooding and gas injection (Martins et al., 1995; Verdon and Stork, 2016) or during carbon geological

storage (Clavijo et al., 2022; Fu et al., 2017). Similarly, thermally-induced fractures are caused due to the injection of cold or hot fluids into the reservoir's rock, creating a thermal shock. The rapid change in temperature causes the rock to fail and develop fractures (Charlez et al., 1996; Martins et al., 1995; Svendson et al., 1991).

During hydraulic fracturing in unconventional reservoirs, proppant placement, with the use of carrying chemicals such as polymers, may cause pore-clogging at the fracture walls, leading to damaging skin and a reduction in matrix-fracture transmissibility (Cramer, 2005; Li et al., 2016; Sugar et al., 2021, 2020). This phenomenon is known as fracture face skin (referred to by $S_f$) (Cinco-Ley, 1981a; Cinco-Ley and Samaniego V, 1977). The fracture face skin phenomenon may also occur in other field applications, such CO$_2$ storage by mineralization (Oelkers et al., 2022, 2021; Omar et al., 2021; Qiao et al., 2021), geothermal

<footnote><sup>a</sup> Physical Science and Engineering Division, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia</footnote>

<footnote><sup>b</sup> Institut Terre et Environnement de Strasbourg, Université de Strasbourg, CNRS, ENGEES, UMR 7063, 67084, Strasbourg, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: hussein.hoteit@kaust.edu.sa (H. Hoteit).</footnote>
<table><thead><tr><th colspan="2">Nomenclature</th></tr></thead><tbody><tr><td>B</td><td>Formation volume factor</td></tr><tr><td>b<sub>s</sub></td><td>Length of the damaged zone</td></tr><tr><td>C<sub>D</sub></td><td>Dimensionless wellbore storage</td></tr><tr><td>c<sub>eff</sub></td><td>Total fracture compressibility</td></tr><tr><td>c<sub>t</sub></td><td>Total compressibility</td></tr><tr><td>F<sub>CD</sub></td><td>Fracture dimensionless conductivity</td></tr><tr><td>h</td><td>Formation thickness</td></tr><tr><td>k<sub>m</sub></td><td>Matrix permeability</td></tr><tr><td>k<sub>f</sub></td><td>Fracture permeability</td></tr><tr><td>k<sub>fi</sub></td><td>Initial fracture permeability</td></tr><tr><td>k<sub>f min</sub></td><td>Minimum fracture permeability</td></tr><tr><td>k<sub>s</sub></td><td>Permeability of the damaged zone</td></tr><tr><td>M</td><td>Number of fracture segments</td></tr><tr><td>N</td><td>Number of fractures</td></tr><tr><td>p</td><td>Formation pressure</td></tr><tr><td>p<sub>D</sub></td><td>Dimensionless pressure</td></tr><tr><td>p<sub>i</sub></td><td>Initial pressure</td></tr><tr><td>p<sub>f</sub></td><td>Fracture pressure</td></tr><tr><td>p&#773;<sub>D</sub></td><td>Formation pressure in Laplace space</td></tr><tr><td>p&#773;<sub>wD</sub></td><td>Pressure at the wellbore in Laplace space</td></tr><tr><td>q</td><td>Flow rate at the subsurface</td></tr><tr><td>q<sub>f</sub></td><td>Flux from the matrix formation to the fracture</td></tr></tbody></table><table><tbody><tr><td>q&#773;<sub>fD</sub></td><td>Flux from the matrix formation to the fracture in Laplace space</td></tr><tr><td>s<sub>D</sub></td><td>Laplace space variable</td></tr><tr><td>S<sub>w</sub></td><td>Wellbore skin factor</td></tr><tr><td>S<sub>f</sub></td><td>Fracture face skin</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>t<sub>D</sub></td><td>Dimensionless time</td></tr><tr><td>w<sub>f</sub></td><td>Fracture width</td></tr><tr><td>x</td><td>x coordinate</td></tr><tr><td>x<sub>f</sub></td><td>Half-length fracture</td></tr><tr><td>y</td><td>y coordinate</td></tr><tr><td>γ</td><td>Formation permeability modulus</td></tr><tr><td>ΔL<sub>fD</sub></td><td>Fracture segment dimensionless length</td></tr><tr><td>μ</td><td>Fluid viscosity</td></tr><tr><td>φ<sub>m</sub></td><td>Matrix porosity</td></tr><tr><td>φ<sub>f</sub></td><td>Fracture porosity</td></tr><tr><td colspan="2"><strong>Subscripts</strong></td></tr><tr><td>D</td><td>Dimensionless</td></tr><tr><td>f</td><td>Fracture</td></tr><tr><td>i</td><td>Initial condition</td></tr><tr><td>m</td><td>Matrix</td></tr><tr><td>min</td><td>Minimum condition</td></tr></tbody></table>

energy extraction from fractured reservoirs that may encounter mineral precipitation (Snæbjörnsdóttir et al., 2020), and waterflooding where scale precipitation at fracture walls may reduce well injectivity (Hajir-ezaie et al., 2017).

Well deliverability in a fractured reservoir is controlled by the fracture network and the matrix-fracture transmissibility. The matrix formation interfacing the fractures may experience a reduction in permeability due to completion or scale precipitation. This reduction is known as fracture face skin, an undesirable condition that creates an additional pressure drop in the system, reducing the productivity or injectivity of the well. The operational conditions and completion practices may be optimized to reduce the fracture face skin. Thus, identifying the type and degree of damage is essential.

During drilling, mud particles may invade the reservoir and induce formation damage. In most formations, excluding highly permeable layers or extreme overbalance pressures, the damaged region is often confined close to the wellbore. However, during hydraulic fracturing, the proppant is placed away from the wellbore to keep the fracture open after the flowback of the fracturing fluid. Embedment into plastic formations of these particles causes the extrusion of the proppant into the formation face, reducing the fracture's flow efficiency. During flowback operations, complete fracturing fluid recovery is rarely achieved, which is caused by fluid retention and loss into the fracture-matrix region (Albattat et al., 2022; Wang and Leung, 2015). In some situations, fluid retention creates a zone of high-water saturation, reducing the reservoir fluid's relative permeability, and causing a decrease in well productivity. This phenomenon can be modeled by the use of equivalent skin, as impairments affecting flow across the fracture walls are accounted for as fracture face skin (Romero et al., 2003).

Various field observations and experimental and modeling studies reported the implications of fracture face skin on flow behavior by many researchers (Adegbola and Boney, 2002; Behr et al., 2006; Gdanski et al., 2005), and experimental work by (Ahmed et al., 1979; Cramer, 2005; Gdanski et al., 2005; Lynn and Nasr-El-Din, 1998; Reinicke et al., 2010; Romero et al., 2003; Roodhart et al., 1988). Fracture face skin often introduces a negative effect on well productivity, where the pressure drop increases with the increase in the skin (Adegbola and Boney, 2002;

Cinco and Samaniego, 1981). Some extreme cases of fracture face skin can cause dramatic degradation of well productivity (Cramer, 2005).

Hydraulic fractures (HF) provide high permeability conduits that enhance well deliverability, which is determined by the fracture hydraulic properties and fracture-matrix crossflow. Several authors reported field experiences indicating degradation of fracture conductivity from its initial level (Clarkson et al., 2013; Ostensen, 1986; Pedrosa, 1986; Zoback and Byerlee, 1975). As the reservoir and fracture pressure deplete during production, the effective stress increases, causing compaction in the fracture system. The resulting degradation in fracture conductivity has been modeled as a function of the pressure drop in the fracture system (Guo et al., 2021; X He et al., 2021; X. He et al., 2020a). On the other hand, the presence of fracture face skin can reduce the fracture aperture and the effective flow area further, leading to an additional pressure drop inside the fracture. Reciprocally, the resulting incremental pressure drop may accelerate the closure process of the fracture. The interconnected effects of fracture closure and fracture face skin have not been considered in the current models in the literature, which is the motivation of this work.

Different techniques are used to estimate fracture properties during and after the hydraulic fracturing procedure, such as micro-seismic, flowback analysis, and borehole imaging (Bennett et al., 2005; Clarkson, 2013; Xupeng He et al., 2021b). One of the most common approaches to assessing fracture flow behavior is by using analytical tools and numerical simulations (Cipolla et al., 2010; Clarkson et al., 2012; Gupta et al., 2012; Lawal et al., 2013; Sun et al., 2015). Numerical simulations require significant computational time and detailed geological information that may not always be available. On the other hand, analytical modeling, based on Pressure Transient Analysis (PTA), is a proven approach adopted in the industry to provide well diagnostics and reservoir parameter estimations with uncertainties (Hao et al., 2015; Hawkes et al., 2013; Xue et al., 2016). Different analytical and semi-analytical PTA models have been developed for fractured wells. Analytical models often consider fractures as parallel planes with a constant effective fracture aperture and simplified fracture geometries (Albattat and Hoteit, 2020; Gringarten et al., 1974; Wong et al., 1986). However, fractures tend to close during production, leading to altered
hydraulic and geometric properties (Hawkes et al., 2013; Xupeng He et al., 2021a). Semi-analytical models with superposition techniques allow modeling fracture closure for unstructured fracture geometries (Jiang et al., 2019a; Yang et al., 2015). To model the effect of fracture closure, the rate of change of fracture conductivity should be incorporated into the fracture governing equations. Previous simulation studies showed that fracture aperture changes significantly at early production time (Liao et al., 2019). Additionally, Chen and Ye (D. Chen et al., 2015) performed loading experiments on shale samples and showed that the effective fracture conductivity is related to the transient fracture aperture under stress conditions, where higher initial compressibility is associated with a larger declining rate of conductivity. Despite the importance of fracture effective compressibility at early times, fractured well analytical and semi-analytical models do not study the effect of fracture compressibility and fracture conductivity reduction. Furthermore, they do not incorporate the impact of fracture face skin.

To model the pressure response for a well with multiple fractures, Yang et al. (2015) proposed a Laplace-space semi-analytical solution based on dividing the fracture system into several blocks and treating each fracture segment as a slab-source term inside the reservoir. In this approach, the hydro-mechanical governing equations in the fracture are discretized similarly to in the Samaniego-Meng solution (Cinco-Ley and Meng, 1988). Later, the permeability modulus property was incorporated to consider pressure-dependent conductivity (Zhang et al., 2014). The fracture equation assumed incompressible flow inside the fracture. The semi-analytical model incorporates the discretized fracture equation and an analytical solution for the reservoir system. The analytical solution is obtained in Laplace space using Newman's principle (Chen et al., 1991; Ozkan and Raghavan, 1991a; Zhang and Yang, 2018). The resulting analytical solution is complex as it consists of an infinite series of hyperbolic functions (Ozkan and Raghavan, 1991b), provided numerical approximations to tackle the numerical divergence issue. With the discretized fracture model and the analytical solution, a system of equations is generated and solved for the pressure response at the wellbore. The system of equations is solved in Laplace space using an iterative method until convergence is achieved for the change in permeability modulus (Jia et al., 2017).

Similarly, Wu et al. (2018) developed a semi-analytical solution for tight-gas fractured wells using pseudo-pressure (Al-Hussainy et al., 1966) to consider the pressure dependency on gas properties and fracture conductivity. In this solution, the fractured system is discretized and coupled with the matrix analytical solution assuming pressure continuity. The matrix solution incorporates the use of a point source solution in Laplace space integrated by the length of each fracture segment (Houze et al., 1988). With the analytical solution and fracture segmentation, a system of equations is generated. The equations are solved for the pressure and the rate for each fracture segment. This model assumes incompressible fluid flow inside the fracture Jia et al. (2017). This fracture discretization is useful for representing complex fracture geometries, including fracture branching, fracture intersections, and non-parallel fractures.

Different approaches have been proposed in the literature to include unstructured fracture geometry to model flow in complex fracture networks. The commonly used approach is to subdivide the fracture network into slabs and assume 1D flow within the slabs (Chen et al., 2018; Z. Chen et al., 2015; Yu et al., 2016; Zhou et al., 2014). Ezulike and Dehghanpour (2014) incorporated the effect of microfractures on wellbore performance using a quadrilinear flow model for natural and hydraulic fracture networks. The solutions are presented in the Laplace Space, dividing the problem into 1D elements at the matrix-fracture interface. Additionally, fluid flow between fracture blocks is accounted for by the use of slab-source functions in the fracture slabs and the superposition principle (Hwang et al., 2013; Lin and Zhu, 2014). The effect of nonplanar fractures, interference, and multiphase flow in semi-analytical models has been incorporated (He et al., 2015; Meng et al., 2020; Yu et al., 2016). Moreover, some models include hydraulic

and natural fracture geometries in the wellbore response (Chen et al., 2018). Clarkson et al. (2013) proposed using a well skin to capture the effect of dynamically-varying fracture conductivity. This approach, however, does not consider the fracture face skin. The effect of partially closed fractures has been considered in various PTA models (Y. He et al., 2020; Huang et al., 2017; Luo et al., 2021, 2020; Qin et al., 2019), where the fracture in segments with lower conductivity for partially propped segments. Soliman (1986) presented a simple mathematical model assuming bilinear flow for a fracture with non-uniform conductivity. Huang et al. (2017) proposed a semi-analytical model based on pseudo pressure for tight gas reservoirs with varying fracture conductivity. These models capture and study the pressure behavior caused by unpropped fracture sections.

Furthermore, models for flowback analysis assuming two-phase flow in complex networks have been proposed (Jia et al., 2018). However, the combined effect of fracture face skin and fracture closure has not yet been addressed, which is the focus of this work. We should note that the fracture face skin should not be confused with the well convergence-skin that has been proposed by many authors to capture the effect of well deliverability degradation because of fracture closure (Bello and Watenbarger, 2010a, 2010b; Clarkson et al., 2013; Zhang and Yang, 2021).

Jiang et al. (2020) incorporated the stress-sensitive matrix permeability, which considers the decrease in fracture conductivity and allows fracture compressibility. However, none of the above-mentioned methods considers the effect of fracture face skin, which is a limitation since it impacts the calculations of the bottom-hole pressure and fracture closure. Wang et al. (2019) presented a method to model fracture aperture during water injection and closure after ceasing injection. This method, however, only considers a single bi-wing fracture with constant matrix-fracture flow along the fracture and assumes a minimum fracture permeability of zero. Other researchers (Wu et al., 2018; Zhang and Yang, 2018) considered wellbore skin but not the actual mechanism of fracture face skin. In this work, we show that the wellbore skin cannot represent the effect of face fracture skin since the skin is not directly incorporated in the boundary condition of the matrix-fracture interface. Moreover, the face fracture skin includes an additional pressure drop to the fracture system, affecting fracture conductivity. The identification, avoidance, and quantification of fracture face skin are crucial in fractured wells, where the fractures provide the primary conduits of flow. However, currently, there is no PTA method that considers complex fracture geometries, fracture closure effects, and the effect of fracture face skin. In this work, we incorporated the fracture face skin into a PTA model that can be used for wells with multiple unstructured fractures. Additionally, the model includes fracture closure and the interaction with fracture face skin.

In this work, we model fractured wells' behavior considering the combined effects of fracture closure and fracture face skin. To generate the pressure response of a fractured well with fracture closure effects and fracture face skin, we divide the reservoir flow regions into three: matrix, skin, and fracture regions. Each region has its respective flow governing equations and is coupled using pressure continuity boundary conditions. The fracture face skin is incorporated as a thin layer of reduced permeability. The fracture's governing equation incorporates the fracture pore volume change and the conductivity reduction due to closure. Superposition in space is used to consider complex fracture geometries. The analytical solution in the matrix is coupled with the fractures by assuming a line-source function for each fracture segment. We use numerical simulations (CMG, 2019) and other analytical methods to verify our model. A field case is also considered to demonstrate the applicability of our model. Using the developed model, we analyze the pressure signals of fracture closure at early, middle, and late times. The interaction of face fracture skin in fracture closure is also discussed.

This paper is organized as follows. We first review the governing equations for compressible single-phase fluid flow in porous media, including a hydraulic fracture. The proposed solution method is then
explained in detail, including the solution method in the matrix, fracture, damage zone, and coupling. Several test cases are then presented to verify and demonstrate the applicability of the proposed model. We end with a discussion about fracture closure detection and conclusions.

## 2. Modeling approach

We subdivide the reservoir domain into three regions corresponding to the matrix, fracture (HF), and the damage zone corresponding to the fracture face skin. Fig. 1 shows a schematic representation of the configuration, where Regions I, II, and III correspond to the HF, matrix, and fracture skin system, respectively. The well is connected to the HF, and the matrix region is connected with the HF through the skin region. The main assumptions of the model are.

- Infinite acting, homogenous and isotropic reservoir with a constant pay thickness.

- The porous medium is saturated with a slightly compressible fluid and constant viscosity.

- The fully penetrated HF's are the pathway to convey fluids from the reservoir to the wellbore.

- The well is operating at a constant production rate.

Using the mass conservation equation and Darcy's law to model single-phase compressible fluid flow in 2D porous media, the governing pressure diffusivity equation is given by (Lee, John, John B. Rollins, 2003):

$$ \frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2} = \frac{\varphi_m \mu c_t}{k_m} \frac{\partial p}{\partial t}, \quad (1) $$

with the initial and boundary conditions,

$$ p(t=0) = p_i, $$

$$ \lim_{x \to \infty} p(t) = p_i, \lim_{y \to \infty} p(t) = p_i, $$

where $\varphi_m$ is the porosity, $\mu$ is the fluid viscosity, $c_t$ is the total matrix compressibility, $k_m$ is the formation permeability, $p$ is the pore pressure, $t$ is time, and $p_i$ is the initial reservoir pressure.

To develop a solution procedure, the governing diffusivity equations are applied separately to the three regions corresponding to the

hydraulic fracture (Region I), the matrix (Region II), and the damage zone (Region III), as illustrated in Fig. 1. The governing equations are then coupled by imposing the continuity of pressure and fluid flow at the region boundaries. The solution procedure is presented in the following subsections, with additional details in the appendices.

### 2.1. Diffusivity model in the matrix

In 2D space, the dimensionless form of Eq. (1) describing the diffusivity equation in the matrix is written as:

$$ \frac{\partial^2 p_{mD}}{\partial x_D^2} + \frac{\partial^2 p_{mD}}{\partial y_D^2} = \frac{\partial p_{mD}}{\partial t_D}, \quad (2) $$

The dimensionless variables are defined by:

$$ p_{mD} = \frac{2\pi k_m h (p_i - p_m)}{q B \mu}, \quad (3) $$

$$ x_D = \frac{x}{x_f}, y_D = \frac{y}{x_f}, t_D = \frac{k_m t}{\varphi_m \mu c_m x_f^2}, \quad (4) $$

where $h$ is the reservoir thickness, $q$ is the surface production rate, $B$ is the formation volume factor describing the fluid unit conversion between reservoir and surface conditions, $x_f$ is the fracture's half-length (see Fig. 1), $p_{mD}$ represents the dimensionless pressure in the matrix (Region I), and $x_D$, $y_D$, and $t_D$ are the dimensionless space and time, respectively. The full derivations of Eq. (1) are provided in Appendix A.

The initial and boundary conditions for an infinite acting reservoir (i. e., in Region II) which contains a fracture with a face skin are written in terms of the dimensionless variables as follows:

$$ p_D(x_D, y_D, t_D = 0) = 0, $$

$$ \lim_{x_D \to \infty} p_{mD}(t_D) = 0, \lim_{y_D \to \infty} p_{mD}(t_D) = 0. $$

At the matrix-skin interface (i.e., Region II-Region III where $y_D = 0$), we impose the continuity of pressure and flow (see Fig. 2), such that,

$$ p_{mD}(y_D = 0) = p_{fD} - q_{fD}S_f, \quad (5) $$

where, $S_f = \frac{\pi b_i}{2x_f} \left( \frac{k_m}{k_i} - 1 \right)$ represents the fracture skin factor, and $q_{fD} = \frac{2q_f x_f}{q_B}$. Note that Eq. (5) incorporates the fracture skin effect and couples

<figure><img src="image_4.png" /><figcaption>Fig. 1. A single hydraulic fracture (HF) illustrated in 3D (left), and in 2D (right) with dimensions and three flow regions corresponding to the matrix, HF, and the damage (skin) zone.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 2. Top 2D view of a horizontal fracture with a fracture face skin (a), and illustration of the pressure profile within the HF, damage zone, and matrix (b), which reflects the continuity of pressure between the HF and the matrix through the damage zone (skin), where a steady-state (linear) pressure is assumed.</figcaption></figure>

the flow between the matrix and the fracture. Additionally, this coupled equation is valid for each fracture segment allowing to consider non-uniform fracture face skin along the fracture. For more details, see Appendix B.

## 2.2. Diffusivity model in the fracture

Consider a HF with a damage zone (Region III) corresponding to the fracture face skin, as illustrated in Fig. 2. The fracture face skin is assumed to be a narrow zone relative to the matrix size with an average thickness $b_s$ and altered permeability $k_s$. We include the effect of fracture face skin into the internal matrix boundary condition. With a relatively thin damage zone compared to the matrix size, we assume steady-state pressure conditions within Region III (see Fig. 2b). Therefore, the extra pressure drop in the damaged zone can be modeled by Eq. (5), which ensures the continuity of the pressure within the skin zone between the matrix and the HF (see Appendix B and C).

Assuming linear flow inside the fracture, the diffusivity equation in the fracture system is given by (see Appendix B):

$$ \frac{\partial^2 p_f}{\partial x^2} - \frac{\pi q_f}{k_f w_f h} = \frac{\varphi_f \mu c_f}{k_f} \frac{\partial p_f}{\partial t}, \qquad (6) $$

with the following initial boundary conditions:

$$ p_f(t=0) = p_i, $$

$$ \left. \frac{\partial p_f}{\partial x} \right|_{x=x_f} = 0, \left. \frac{\partial p_f}{\partial x} \right|_{x=x_w} = - \frac{q_B \mu}{k_f w_f h}, $$

where $\varphi_f$ is the fracture porosity, $x$ is the distance along the fracture, $c_f$ is the total fracture compressibility, including the change in fracture aperture, $w$ is the fracture width, $p_f$ is the fracture pore pressure, $k_f$ is the fracture permeability, which is a function of pressure (see Eq. (24)), $q_f$ is the rate coming from the matrix to the fracture per unit length, and $p_i$ is the initial pressure. The effective fracture compressibility incorporates the change in fracture aperture due to closure in addition to fluid and

rock compressibility, i.e., $c_{ft} = \frac{1}{w} \frac{\partial w}{\partial p_f} + c_f$. The fracture diffusivity equation can be transformed to dimensionless form with the following variables (see Appendix B):

$$ p_{fD} = \frac{2\pi k_m h (p_i - p_f)}{q_B \mu}, \qquad (7) $$

$$ F_{CD} = \frac{k_f w}{k_m x_f}, F_{CD}^* = F_{CD} \delta_D, \qquad (8) $$

$$ q_{fD} = \frac{2q_f x_f}{q_B}, \gamma_D = \frac{q_B \mu \gamma}{2\pi k h}, \qquad (9) $$

$$ \delta_D = \frac{F_{CDmin}}{F_{CDi}} + \left(1 - \frac{F_{CDmin}}{F_{CDi}}\right) e^{(-\tau_D p_{fD})}, \qquad (10) $$

where $\gamma$ is the permeability modulus of the fracture, $\delta_D$ is used to model the effect of conductivity reduction (Zhang et al., 2014). The parameter is obtained using the definition of dimensionless variables. It assumes a threshold in the permeability reduction and incorporates the permeability modulus (Berumen and Tiab, 1997). Further details can be found in the section on pressure-dependent fractures. The dimensionless form of Eq. (6) is written as:

$$ \frac{\partial^2 p_{fD}}{\partial x_D^2} - \frac{\pi q_{fD}}{F_{CD}} = \frac{1}{\eta_D} \frac{\partial p_{fD}}{\partial t_D}, \qquad (11) $$

with the initial and boundary conditions,

$$ p_{fD}(x_D, t_D = 0) = 0, \frac{\partial p_{fD}}{\partial x_{D,x_D=1}} = 0, $$

$$ \frac{\partial p_{fD}}{\partial x_{D,x_D=0}} = - \frac{\pi}{F_{CD}}. $$

To solve the pressure within the HF, we use a semi-analytical approach where the HF is discretized into 2M subblocks, as shown in Fig. 3.

## 2.3. Matrix-fracture coupling

A semi-analytical method is used to couple and solve the diffusivity equations in the matrix and fracture regions. The semi-analytical solution in the matrix is solved within Laplace space, and the results are converted into time-space, while the solution in the fracture is performed numerically. In Laplace space, the diffusivity equation with boundary conditions in the matrix is given by:

$$ \frac{\partial^2 \bar{p}_D}{\partial x_D^2} + \frac{\partial^2 \bar{p}_D}{\partial y_D^2} = s_D \bar{p}_D, \qquad (12) $$

$$ \lim_{x_D \to \infty} \bar{p}_D(s_D) = 0, \lim_{y_D \to \infty} \bar{p}_D(s_D) = 0, $$

$$ \bar{p}_D(y_D = 0) = \bar{p}_{fD} - \bar{q}_{fD} S_f. \qquad (13) $$

Where $s_D$ is the Laplace transformation variable. The matrix diffusivity equation can be obtained by considering each HF subblock as a source/ sink term (see Fig. 3), following the solution method proposed by (Jiang et al., 2021). Assuming N fractures divided each into 2M fracture segments each of length $\Delta L_{fD}$, the pressure solution in Laplace space corresponds to each fracture block, labeled as $(i,j)$, is given by:

$$ \bar{p}_{Di,j} = \frac{1}{2} \sum_{j=1}^{2NM} \bar{q}_{fjDj} \int_{-\Delta L_{fD}/2}^{\Delta L_{fD}/2} K_0 \left[ \sqrt{s_D} \sqrt{(x_D - x_{Di} - u)^2 + (y_D - y_{Di})^2} \right] du \qquad (14) $$

Equation (14) calculates the dimensionless pressure response at each
<figure><img src="image_6.png" /><figcaption>Fig. 3. A 2D top view of the HF with the surrounding face skin and matrix zones, where the HF is discretized into 2M subblocks; the arrows indicate the fluid flow directions.</figcaption></figure>

fracture segment (x<sub>D</sub>,y<sub>D</sub>), where the superposition in space is implemented to account for the interference of the different fracture blocks located at (x<sub>wDib</sub>, y<sub>wDib</sub>). With this superposition principle, the pressure solution can be obtained for multiple fractures with non-uniform configurations (see Fig. 4). The pressure response is obtained by solving a system of equations, where the matrix and fracture solutions can be arranged in a matrix form. The pressure response at spatial location (x<sub>Di</sub>, y<sub>Di</sub>) can be obtained as follows:

$$A\bar{q}_{\text{FD}} - I\bar{p}_{\text{D}} = \bar{0}, \qquad (15)$$

where A = [A<sub>ij</sub>] is an invertible (N × 2M, N × 2M) matrix consists of submatrices A<sub>ij</sub> for each fracture segment i and using the superposition principle, and q̄<sub>FD</sub> and p̄<sub>D</sub> correspond to the flux and pressure variables, given by:

$$A_{i,j} = \frac{1}{2} \int_{-\Delta L_{\text{FD}}/2}^{\Delta L_{\text{FD}}/2} K_0 \left[ \sqrt{s_D} \sqrt{(x_D - x_{wDi} - u)^2 + (y_D - y_{wDi})^2} \right] du + S_f, \quad (16)$$

$$\bar{p}_{\text{D}} = [ \bar{p}_{D,1} \quad \bar{p}_{D,2} \quad \dots \quad \bar{p}_{D,N \times 2M} ]^{\text{T}},$$

$$\bar{q}_{\text{FD}} = [ \bar{q}_{fD,1} \quad \bar{q}_{fD,2} \quad \dots \quad \bar{q}_{fD,N \times 2M} ]^{\text{T}},$$

In Eq. (15), I is the identity matrix with dimensions (N × 2M, N × 2M), 0̄ is the vector zero with dimension N × 2M. Note that Eq. (16) incorporates the effect of the fracture face skin by imposing the boundary condition defined in Eq. (13) for each fracture segment.

Note that both variables p̄<sub>D</sub> and q̄<sub>FD</sub> are unknowns in Eq. (15), and therefore another system corresponding to the governing equations in the fractures is needed to close the solution procedure.

The governing equations in the fractures using Laplace space can be derived from the dimensionless diffusivity equation given in Eq. (11), such that,

$$\frac{\partial^2 \bar{p}_{fD}}{\partial x_D^2} - \frac{\pi \bar{q}_{fD}}{F_{CD}^*} = \frac{s_D}{\eta_D} \bar{p}_{fD}, \qquad (17)$$

<figure><img src="image_7.png" /><figcaption>Fig. 4. Nonuniform configuration of multiple fractures (N), where each fracture is subdivided into a number of segments (2 M), resulting in N × 2 M fracture blocks.</figcaption></figure>

$$\frac{\partial \bar{p}_{fD}}{\partial x_{D, x_D=0}} = - \frac{\pi}{F_{CD}^* s_D} \frac{\partial \bar{p}_{fD}}{\partial x_{D, x_D=1}} = 0, \qquad (18)$$

$$\eta_D = \frac{k_f \varphi_m c_{om}}{k_m \varphi_f c_{if}}. \qquad (19)$$

As previously discussed, we use a finite difference discretization for the hydraulic fractures, where each fracture is subdivided into fracture blocks (see Fig. 3). The material balance equation is solved by imposing the continuity of the fluxes, in particular, the matrix-fracture flux, which couples the flow between the HF and formation matrix systems, including the fracture face skin.

The resulting (N × 2M, N × 2M) linearized system corresponding to the flow equations in the HF becomes,

$$F\bar{q}_{\text{FD}} + E\bar{p}_{\text{FD}} + B\bar{p}_{\text{wD}} = \bar{0}, \qquad (20)$$

with,

$$\bar{p}_{\text{FD}} = [ \bar{p}_{fD,1} \quad \bar{p}_{fD,2} \quad \dots \quad \bar{p}_{fD,N \times 2M} ]^{\text{T}}, \qquad (21)$$

where 0̄ is defined as in Eq. (15), F is a diagonal matrix with dimensions (N × 2M, N × 2M), and diagonal values correspond to ω<sub>Di</sub>, E is a diagonal matrix with dimensions (N × 2M, N × 2M), and B is a vector of size N × M. A detailed description of the procedure is given in Appendix D.

## 2.4. Solution procedure

The solution procedure consists of coupling the diffusivity equations that account for the fracture closure, including pressure-dependent permeability (see Eq. (24)), fracture aperture reduction, and fracture face skin effect described by permeability alteration around the fracture walls. The resulting systems are coupled by the continuity of pressure and mass flux between the matrix and fracture regions. The fracture face skin is accounted for as a boundary condition for the matrix system. The resulting linear system of equations from the matrix and fracture regions are provided in Eq. (15) and Eq. (20), respectively. By integrating the two systems, the final linear system of equations becomes,

$$\begin{bmatrix} A & -I & 0 \\ F & E & B \\ 0 & B' & C \end{bmatrix} \begin{bmatrix} \bar{q}_{fD} \\ \bar{p}_{fD} \\ \bar{p}_{wD} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ D \end{bmatrix} \qquad (22)$$

The pressure variables in the Laplace space are calculated by solving the above system of equations. Then, the solution in the time-space is calculated numerically using the Stehfest algorithm (Stehfest, 1970). A detailed description of the solution procedure and the structure of the linear system in Eq. (22) are provided in Appendix D.

## 2.5. Pressure-dependent fractures

The mechanism of pressure-dependent fracture conductivity in the fractures is incorporated using a dimensionless formation permeability modulus (γ), which describes the change in permeability with respect to the change of pressure (Berumen and Tiab, 1997; X He et al., 2021; X. He et al., 2020b; Pedrosa, 1986). The permeability modulus gives an
exponential reduction of permeability with pressure decline. The pressure-dependent fracture permeability ($k_f$) at a given fracture pressure ($p_f$) is given by:

$$k_f(p_f) = k_{fi}e^{-\gamma(p_i-p_f)}, \qquad (23)$$

where $k_{fi}$ is the initial fracture permeability and $p_i$ is the initial (or reference) pore pressure in the fracture. To consider a threshold in permeability reduction, we adopt the model proposed by (Zhang et al., 2014), where the pressure-dependent permeability is modified to,

$$k_f(p_f) = k_{f \min} + (k_{fi} - k_{f \min})e^{-\gamma(p_i-p_f)}, \qquad (24)$$

where $k_{f \min}$ is a lower limit of fracture permeability. Note that Eq. (24) reduces to Eq. (23) when $k_{f \min} = 0$.

Using Eq. (24) and the definition of the dimensionless variables, the relationship given by Eq. (10) is obtained. The maximum and minimum permeability values can be obtained from experimental tests (Aguilera, 1999; Jiang et al., 2019b; Pedrosa, 1986). These laboratory tests account for the proppant crushing an embedment, and permeability-damage mechanisms are performed to explain the loss of productivity (Palisch et al., 2007). The obtained permeability values are used along with the pressure transient data for identifying the conductivity reduction.

## 3. Model verification

### 3.1. Test case 1: fracture face skin effect with no fracture closure

We compare our proposed model with a simplified analytical model by Cinco-Ley and Samaniego (Cinco-Ley and Samaniego V, 1977) to capture the fracture skin effect. The model assumes incompressible fluid flow inside a fracture with a bi-wing configuration and neglects fracture closure. The Cinco-Ley-Samaniego solution considers a vertical well with a vertical fracture with a constant thickness and a surrounding zone of diminished permeability. The Cinco-Ley-Samaniego's problem is a particular case for our model. Fig. 5 shows the results for the dimensionless pressure ($p_D$) versus dimensionless time ($t_D$) for the two solutions, which reflect excellent consistency. This problem shows that as the fracture face skin ($S_f$) increases an additional pressure drop is generated from early to late times. As a result, it is expected that the incorporation of fracture face skin can affect the fracture closure behavior since the closure is modeled as a function of the pressure drop generated. This scenario is investigated in the following test case.

<figure><img src="image_8.png" /><figcaption>Fig. 5. Comparison of solutions obtained by the proposed model and the Cinco-Ley-Samaniego model, reflecting an increase in dimensionless pressure, $p_D$ (or pressure drop) as the fracture face skin ($S_f$) increases; the results show an excellent match between the two models (Test case 1).</figcaption></figure>

### 3.2. Test case 2: fracture face skin effect with fracture closure

Because of the absence of analytical or semi-analytical models to capture the combined effect of fracture skin and closure, we used numerical simulation to verify our model. A simulation model was built using CMG-IMEX (CMG, 2019) reservoir simulator to verify the effects of fracture closure and fracture face skin. The simulation model includes a single hydraulic fracture. The fluid and rock properties in the matrix and fracture are listed in Table 1. Fig. 6 shows the grid used to perform the simulation where local grid refinement was conducted in the near-wellbore region and around the HF. The model considers a bi-wing fractured well producing at a constant rate. A thin layer of reduced permeability was introduced around the fracture to mimic the skin effect. The fracture face skin is assumed to be constant along the fracture. However, the proposed model allows to incorporate varying fracture face skin along the fracture. Table 2 shows the values of the reduced permeability and the length of the skin zone. The used effective fracture compressibility is adopted from Aguilera (1999) for propped fractures, which is based on experimental data (Jones, 1975; Tkhostov et al., 1970).

Fig. 7 shows the dimensionless pressure response and its derivative obtained from the simulations and the proposed semi-analytical model, which are in excellent agreement. The pressure derivative curve is shown as $p_D/10$ for visualization purposes. Results show that at early times, linear flow is observed with a slope of 1/2 in the derivative curve, followed by a transition period to the radial flow with a 0 slope in the derivative. The incorporation of fracture face skin creates an additional pressure drop at the fracture-matrix interface, as expected. As a result, the pressure derivative shows an increment and a deviation from the formation linear flow. Similarly, the face fracture skin increases the linear flow duration, which is consistent with the results reported in the literature (Valdes-Perez et al., 2013; Wong et al., 1986). In this test case, the pressure-dependent permeability is not included.

### 3.3. Test case 3: fracture closure, fracture face with pressure-dependent permeability

In this example, we use reservoir simulations to verify the proposed model, including the effect of fracture closure, permeability-dependent pressure, and fracture face skin. We use the reservoir parameters provided in Table 1. The face skin is represented by a thin layer of reduced permeability, and the change in fracture permeability is set as a function of pressure. A pressure-permeability look-up table is used to describe the pressure-dependent permeability of the fracture in the CMG simulation model. Fig. 8 shows the combined effect of face fracture skin and fracture closure using simulations and the proposed model. The pressure derivative curve is shown as $p_D/10$ for visualization purposes. For reference, we also show the base-case solution with no fracture closure.

<table><caption>**Table 1**

Parameters used for model verification with the numerical model.</caption><thead><tr><th colspan="2">Reservoir properties</th></tr></thead><tbody><tr><td>matrix permeability, k[md]</td><td>1.0</td></tr><tr><td>matrix porosity, φ[−]</td><td>0.12</td></tr><tr><td>matrix compressibility, c<sub>m</sub>[psi<sup>-1</sup>]</td><td>0.0</td></tr><tr><td>net pay thickness, h[ft]</td><td>82</td></tr><tr><td>initial reservoir pressure, p<sub>i</sub>[psi]</td><td>8770</td></tr><tr><td colspan="2"><strong>Oil properties</strong></td></tr><tr><td>formation Volume Factor, B<sub>o</sub>[rbl/stb]</td><td>1.03</td></tr><tr><td>oil viscosity, μ<sub>o</sub>[cp]</td><td>0.65</td></tr><tr><td>compressibility, c<sub>o</sub>[psi<sup>-1</sup>]</td><td>5.2 × 10<sup>-6</sup></td></tr><tr><td colspan="2"><strong>Hydraulic fracture properties</strong></td></tr><tr><td>fracture half-length, x<sub>f</sub>, [ft]</td><td>500</td></tr><tr><td>initial fracture conductivity, F<sub>CDI</sub>[−]</td><td>20</td></tr><tr><td>dimensionless permeability modulus, γ<sub>D</sub>[−]</td><td>1</td></tr><tr><td>minimum fracture conductivity, F<sub>CDmin</sub>[−]</td><td>0.2</td></tr><tr><td>effective fracture compressibility, c<sub>eff</sub>[psi<sup>-1</sup>]</td><td>3.2 × 10<sup>-4</sup></td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Simulation model showing a single fracture with local grid refinement and a face skin zone around the fracture region (Test case 2).</figcaption></figure>

<table><caption>**Table 2**
Permeability and extension of the damaged (skin) zone surrounding the fracture.</caption><thead><tr><th>Fracture face skin, S<sub>f</sub></th><th>Permeability of skin zone, k<sub>s</sub>[md]</th><th>Extension of the skin zone, b<sub>s</sub>[ft]</th></tr></thead><tbody><tr><td>0.00</td><td>1.00</td><td>—</td></tr><tr><td>0.04</td><td>0.08970</td><td>1.256</td></tr><tr><td>0.10</td><td>0.03794</td><td>1.256</td></tr><tr><td>0.26</td><td>0.01494</td><td>1.256</td></tr><tr><td>0.50</td><td>0.00782</td><td>1.256</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Comparison of solutions obtained by the proposed model and simulations (CMG-IMEX), reflecting an increase in dimensionless pressure, $p_D$ and its derivative, $p'_D$ as the fracture face skin ($S_f$) increases; results show excellent agreement between the two models (Test case 2).</figcaption></figure>

and no face skin. The results show a good agreement between the simulation and the proposed semi-analytical solution. The effect of face fracture skin is mainly observed in the early and middle times, delaying the beginning of the transition flow. The face fracture skin also creates an additional pressure drop, observed as an increment in the dimensionless pressure $p_D$. As a consequence of this skin-induced pressure drop, fracture closure is triggered sooner. This test case shows that the proposed model is able to capture the combined effect of fracture closure and fracture face skin.

### 3.3.4. Test case 4: fracture closure and face skin for a multistage fractured well

To verify the applicability of the proposed model for a multistage hydraulic fracture well, a simulation model was built for a 5-stage hydraulic fracture well (see Fig. 9). The properties of the model are the

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Comparison of dimensionless pressure curves obtained by simulations and the proposed model, where the combined effect of fracture closure and fracture face skin are considered (Test case 3).</figcaption></figure>

same for the fractures and the matrix, as listed in Table 1. The fracture spacing is set at 1000 ft. Each fracture is surrounded by a thin skin layer of reduced permeability corresponding to a face skin of, $S_f = 0.30$.

Fig. 10 compares the numerical solution and the proposed semi-analytical solution for the multistage fracture cases, which shows excellent agreement between the two models. For reference, we also show the solution without including the fracture face skin. The effect of the fracture face skin is more evident during the linear flow and the bilinear flow. These two flow regimes are the most prevalent in hydraulically fractured wells, where low reservoir permeability requires long production times to develop pseudo-radial or radial flow regimes. The fracture face skin induces an additional pressure drop, resulting in increasing the slope of the dimensionless pressure and its derivative. At later times, the pressure response converges in the pseudo-radial behavior. In this example, the effect of the fracture face skin is not noticeable at the late time. Further details are provided in the following Discussion session.

To verify the applicability of the proposed model for a horizontal multistage hydraulic fractured well, a simulation model was built for an 11-stage hydraulic fracture well in a low permeability reservoir. The properties of the model are listed in Table 3. The fracture spacing is set at 50 ft. Each fracture is surrounded by a thin skin layer of reduced permeability corresponding to a face skin of $S_f = 0.30$.

Fig. 11 compares the numerical solution and the proposed semi-analytical solution for the multistage fracture cases, which shows excellent agreement between the two models. For reference, we also show the solution without including the fracture face skin. The effect of the fracture face skin is more pronounced during the linear flow and the bilinear flow. In this example, the effect of the tip of the fracture is considered (Wong et al., 1986). Since the fracture face skin impedes the interaction between the fracture and the matrix system, the bilinear flow is delayed, and the pressure response resembles a closed system since no-flow boundaries are imposed on the fractures' tips. In this case, the effect of fracture face skin is more evident as linear, and bilinear flow regimes are the most prevalent in hydraulically fractured wells due to the low reservoir permeability. Additionally, a similar behavior in pressure is observed where the fracture face skin induces an additional pressure drop, resulting in an increase in the slope of the dimensionless pressure and its derivative. Following this period, bilinear flow is observed, and at late times, the pressure response shows the effect of the fracture interacting with each other.
<figure><img src="image_12.png" /><figcaption>Fig. 9. Simulation model showing a multistage fractured with local grid refinement and a face skin zone around the fracture region (Test case 4).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Comparison of dimensionless pressure considering the combined effect of fracture closure and fracture face skin for a multistage fractured well using the proposed model and those from numerical simulations (Test case 4).</figcaption></figure>

<table><caption>**Table 3**

Parameters used for model verification with the numerical model.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td><strong>Reservoir properties</strong></td><td></td></tr><tr><td>matrix permeability, k[md]</td><td>0.001</td></tr><tr><td>matrix porosity, φ[ - ]</td><td>0.12</td></tr><tr><td>matrix compressibility, c<sub>m</sub>[psi<sup>-1</sup>]</td><td>0.0</td></tr><tr><td>net pay thickness, h[ft]</td><td>82</td></tr><tr><td>initial reservoir pressure, p<sub>i</sub>[psi]</td><td>8770</td></tr><tr><td><strong>Oil properties</strong></td><td></td></tr><tr><td>formation Volume Factor, B<sub>o</sub>[rbl/stb]</td><td>1.26</td></tr><tr><td>oil viscosity, μ<sub>o</sub>[cp]</td><td>0.65</td></tr><tr><td>compressibility, c<sub>o</sub>[psi<sup>-1</sup>]</td><td>2.1 × 10<sup>-7</sup></td></tr><tr><td><strong>Hydraulic fracture properties</strong></td><td></td></tr><tr><td>fracture half-length, x<sub>f</sub>[ft]</td><td>500</td></tr><tr><td>initial fracture conductivity, F<sub>CDI</sub>[ - ]</td><td>20,000</td></tr><tr><td>dimensionless permeability modulus, γ<sub>D</sub>[ - ]</td><td>1</td></tr><tr><td>minimum fracture conductivity, F<sub>CDmin</sub>[ - ]</td><td>0.2</td></tr><tr><td>effective fracture compressibility, c<sub>eff</sub>[psi<sup>-1</sup>]</td><td>3.2 × 10<sup>-4</sup></td></tr></tbody></table>

## 4. Results and discussion

This section analyzes pressure signals for early, middle, and late times considering the effects of fracture closure (i.e., fracture-effective compressibility and pressure-dependent permeability). Depending on the pressure data collected in the field during the transient pressure test, three-time regimes can be interpreted to estimate if the fracture is exhibiting closure and the closure rate. In general, pressure response at early times is sensitive to effective fracture compressibility, while in middle times, the pressure response is mainly affected by the fracture

<figure><img src="image_14.png" /><figcaption>Fig. 11. Comparison of dimensionless pressure considering the combined effect of fracture closure and fracture face skin for a multistage fractured well using the proposed model and those from numerical simulations in a low permeability reservoir.</figcaption></figure>

conductivity ($F_{CD}$) reduction. For low dimensionless permeability modulus $γ_D$, fracture closure is observed at late times using a semi-log approach.

### 4.1. Early times response

The total effective fracture compressibility, $c_{ft}$, has a significant effect at early times. It covers bilinear flow, and a long linear flow is observed. This behavior was reported in previous works (Berumen and Tiab, 1997; Cinco-Ley, 1981b; Liao et al., 2019; Tiab et al., 1999). (Cinco-Ley, 1981b) described this effect as a fracture with high storage capacity (see Fig. 12). It is important to point out that this effect is unrelated to the fracture storage period since the pressure derivative slope does not correspond to fracture storage. The behavior of a high storage capacity fracture is related to the total compressibility of the fracture system. For a high storage capacity fracture, the derivative slope is $m = 1/2$; meanwhile, the fracture storage slope is $m = 1$ (Wong et al., 1986).

Using the parameters in Table 1, we incorporated the effective fracture compressibility to mimic the behavior of a high storage capacity. Using this configuration, the values obtained for fracture conductivity and dimensionless diffusivity are $F_{CDI} = 20$, $η_D = 500$, corresponding to a fracture with high compressibility.

Fig. 13 shows the comparison of the model using the effective fracture compressibility and neglecting the compressibility of the fracture with $F_{CDmin} = 0.2$. The effect of fracture effective compressibility is visible at early times, and it covers the bilinear flow with a long linear flow and a lower pressure drop. Fracture compressibility does not affect the pressure response at middle and late times. This behavior is
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Typical behavior of dimensionless pressure versus dimensionless time in vertical wells with a hydraulic fracture.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Early-time response with and without including the effective fracture compressibility.</figcaption></figure>

consistent with the results observed for flowback simulation results and field data (Liao et al., 2019; Tiab et al., 1999).

The effect of fracture compressibility becomes more pronounced as the pore volume of fracture formation increases (see Fig. 13). In massive hydraulic fracturing, the fractures' pore volume is high compared to a single vertically fractured well. As a result, fracture closure affects the flowback and initial production period. The early flow period in fractured wells is used as a post-fracturing analysis to estimate fracture properties. Furthermore, initial fracture compressibility can vary according to the proppant placement, mineralogy, and anisotropy (D. Chen et al., 2015). With the incorporation of fracture-effective compressibility, the model can be used to analyze the flowback period and early production time. Additionally, the early time behavior can be used as an indication of the degree of compressibility of the fractured system.

## 4.2. Middle times

Middle times correspond to the late bilinear flow and the beginning of the radial flow. This flow period is mainly dominated by the change in fracture conductivity, but it is barely affected by fracture compressibility (see Fig. 14). The effect of pressure-dependent fracture conductivity is incorporated in the dimensionless permeability modulus ($\gamma_D$). For a high

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Comparison of dimensionless pressure curves from this model and those from numerical simulations considering the effect of minimum fracture conductivity.</figcaption></figure>

value of permeability modulus, the pressure response during the transition flow is affected, generating an increment in the dimensionless pressure and its derivative. If the fracture closes slowly (i.e., low $\gamma_D$), the effect will be reflected in the late-times region. The occurrence of fracture permeability reduction generates an additional pressure drop and generates a 'hump' in the pressure derivative. Using the pressure response and its derivative, the match of pressure data can be achieved (Jiang et al., 2019a). Nevertheless, when the fracture conductivity changes are slow, or the initial fracture conductivity is high, the closure effects are evident in the radial flow. However, they may not be identified during the transition flow.

At middle times, fracture closure affects the conductivity of the fracture. Conductivity reduction results in an increment of the pressure drop in the system. The dimensionless pressure increases when conductivity starts reducing, generating an increment in the dimensionless pressure derivative (see Fig. 14). Once the conductivity reaches the minimum conductivity $F_{CDmin}$, the effect on the derivative becomes less evident. Nevertheless, other researchers (Jiang et al., 2019a) reported that in high-conductivity fractures, no significant stress-sensitivity effects could be found during the transient pressure behavior. However, our results indicate that fracture closure of high-conductivity fractures is visible depending on the minimum fracture conductivity $F_{CDmin}$. Fig. 14 shows the pressure response from the proposed model and simulations considering different values for minimum fracture conductivity $F_{CDmin}$, with $\gamma_D = 3$ for a vertical fractured well. Results show that for 90% loss of fracture conductivity, fracture closure is barely detected. However, as the minimum fracture conductivity decreases, the deviation of the pressure response becomes evident. In other words, fracture closure can be identified in high-conductivity fractures, and it depends on the degree of fracture conductivity reduction and not on the initial fracture conductivity. This behavior becomes more evident in high-conductivity fractures where $F_{CDmin}$ is the indicator to identify fracture closure, regardless of the rate of conductivity reduction. This point will be extended in the following discussion for the late-times response.

## 4.3. Late times

In this section, the effects of fracture closure after the transition period are discussed. When the fracture conductivity changes are slow, or the initial fracture conductivity is high, the effect of fracture closure is evident in the radial flow, which may not be identified during the transition flow. By assuming a high fracture conductivity or a fracture closing at slow rates, the pressure response at early and middle times will not be detected. The effect of pressure-dependent fracture conduc-
tivity is incorporated into the dimensionless permeability modulus ($\gamma_D$). For a high value of permeability modulus and low conductivity fractures, the pressure response during the transition flow is detected by generating an increment in the dimensionless pressure and its derivative. If the fracture closes at a slow rate, the effect will be visible in the late-time region.

To identify fracture closure in high conductivity fractures or at slow closure rates (i.e., low permeability modulus), we propose using a semi-log plot to analyze the pressure response during radial flow. Fig. 15 shows the pressure response for a vertical hydraulic fractured well with initial fracture conductivity, $F_{CDi} = 300$, decreasing its conductivity at different rates. At late times, the pressure response tends to drop to the minimum fracture conductivity $F_{CDmin} = 1$. The deviation in the initial pressure behavior corresponds to the conductivity reduction due to closure. This deviation becomes more evident in a semi-log scale than in a log-log plot.

To demonstrate the use of the semi-log approach, we built a CMG model for a fractured well with equivalent properties of $F_{CDi} = 300$ that reduces its conductivity $F_{CDmin} = 1$ using $\gamma_D = 1.0$. Fig. 15 shows the simulation results and the semi-analytical response for different rates on a log-log scale. As the fracture closure mainly occurs in the radial flow, $p_D$ and $p_D$ in the log-log plot make the fracture's closure rate identification uncertain. The semi-log plot is recommended for the identification of fracture closure (see Fig. 15) since the pressure drop change is more evident in the dimensionless pressure behavior than in its derivative. By analyzing the pressure curves from different fracture conductivities in radial flow, this approach performs better when $F_{CDmin} < 50$. This can be explained considering that in high conductivity fractures

(infinite conductivity fractures), the pressure drop along the HF can be neglected. Once the conductivity is reduced to finite conductivity values, the pressure drop inside the fractures becomes evident. This behavior is related to the concept of effective wellbore radius for a fractured well.

The semi-log approach can be simplified using the concept of effective wellbore radius (Prats et al., 1962). For a hydraulic fractured well with a finite conductivity fracture, the radial flow can be analyzed using the approach presented by Cinco-Ley and Samaniego (Cinco-Ley, 1981b). At late times the fracture behaves like a well with an equivalent radius. Fig. 16 shows the relation between fracture conductivity and effective wellbore radius incorporating a negative skin value $r'_w = r_{weff}^{-s}$. When radial flow is achieved in the well, the dimensionless pressure can be computed using the expression:

$$p_D = \frac{1}{2} \ln \left( \frac{4t_D}{\gamma} \right) + \ln \left( \frac{2r'_w}{x_f} \right) \quad (25)$$

By using this simple expression, radial flow can be analyzed to identify the fracture closure. Using the previous example where a fracture is reducing its conductivity from $F_{CDi} = 300$ to a minimum fracture conductivity $F_{CD} = 1$. Fig. 16 shows the equivalent radius values (red dots) for the simulation model results for initial, final, and current fracture conductivity.

The current fracture conductivity can be estimated using the effective radius and the semi-log approach in radial flow. Fig. 16 shows the pressure response from the simulations, the proposed semi-analytical model, and the effective wellbore radius solution. Using the semi-log approach, the concept of effective radius $r'_w$ is helpful to observe

<figure><img src="image_18.png" /><figcaption>Fig. 15. Comparison of log-log (left) and semi-log (right) approach in radial flow to identify fracture closure.</figcaption></figure>

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 16. Effective wellbore radius implementation to identify fracture closure and used values of effective wellbore radius.</figcaption></figure>

changes in fracture conductivity in time with a less expensive compu-
tational method. Fig. 16 shows that the effective wellbore radius
matches the pressure response computed using the proposed semi-
analytical solution at constant F<sub>CD</sub> = 1 and F<sub>CD</sub> = 300. When closure
becomes evident, the pressure response deviates from the initial effec-
tive wellbore radius and intersects another value of an equivalent radius,
which corresponds to the current fracture conductivity. Fig. 16 also
shows a threshold value of F<sub>CD</sub> = 50 where the effective wellbore radius
becomes practically constant and the dimensionless pressure behavior
becomes almost identical, and for which fracture closure will not be
identified. It is important to note that the effective wellbore radius
considers a constant conductivity along the fracture and the semi-
analytical model considers the change of fracture conductivity in each
fracture segment. Thus, the conductivity obtained using the effective
wellbore radius is an equivalent fracture conductivity for all fracture
segments. One disadvantage of the pseudo-radial approach is the lack of
pressure derivative matching since the pseudo-radial approach deriva-
tive will always be 1/2. It is important to consider that in low permeability
reservoirs, achieving pseudo-radial flow requires long production times.

4.4. *Fracture face skin and closure effects*

As fracture closure is a function of the internal pressure of the frac-
ture system, fracture face skin directly impacts the fracture closure.
Fig. 17 compares the proposed semi-analytical solution for fracture

<figure><img src="image_21.png" /><figcaption>Fig. 17. Dimensionless pressure behavior with and without considering fracture face skin.</figcaption></figure>

<figure><img src="image_22.png" /></figure>

closure, with and without considering face fracture skin for a vertical
fractured well. We consider F<sub>CDi</sub> = 20, F<sub>CDmin</sub> = 0.2, η<sub>D</sub> = 500. The first
increase in pressure drop is related to the face fracture skin, followed by
the pressure drop generated by the fracture closure. As the fracture face
skin increases, the fracture conductivity reduction is triggered earlier, as
appears in Fig. 17.

The presence of a damaged zone surrounding the fracture creates an additional pressure drop, which increases the fracture’s effective stress, triggering the fracture closure. Fig. 18 shows the effect of fracture face skin on the reduction of fracture conductivity. As the value of fracture face skin increases, the conductivity reduces further.

4.5. *Fracture face skin versus wellbore skin*

The current approaches in the literature that use wellbore skin or convergence-skin to capture the effect of fracture face skin are not accurate (Bello and Wattenbarger, 2010a, 2010b; Clarkson et al., 2013; Zhang and Yang, 2021). Wellbore skin and storage can be added using the classical Laplace equation, given by (Van Everdingen and Hurst, 1949). The effect of the wellbore skin and the fracture face skin have different behaviors. The wellbore skin adds an additional pressure drop based on the total well flow rate. On the other hand, the fracture face skin considers the pressure drop in each of the fracture segments that may affect the fracture, as discussed previously. In previous models, once the pressure at the wellbore (p<sub>wD</sub>) is obtained, the skin is added as

<figure><img src="image_23.png" /><figcaption>**Fig. 18.** Fracture conductivity reduction resulting from different fracture face skin.</figcaption></figure>
an additional pressure drop using the original definition of wellbore ski, given by:

$$ \bar{p}_{wD} = \frac{s\bar{p}_{wD} + S_w}{s + C_D s^2 (\bar{p}_{wD} + S_w)}, \quad (26) $$

where $C_D$ is the dimensionless wellbore storage and $S_w$ is the wellbore skin factor. However, this approach does not represent the permeability changes near the fracture region. Fig. 19 shows the results from simulations and the comparison of considering wellbore skin $S_w$ (previous models) and the proposed fracture face skin ($S_f$) for a vertical fractured well. Including face fracture skin in the solution induces an extra pressure drop and generates a 'hump' in the pressure derivative. However, previous models only modify the pressure response at early times, and the resulting pressure derivative does not capture the correct behavior. Moreover, the usual skin approach does not affect the fracture conductivity change, which may be erroneous, considering the conductivity change is a function of the pressure drop inside the fractured system.

# 5. Field example

## 5.1. Water injection well

In oilfields undergoing waterflooding, water injectors may induce formation fracturing due to uncontrol injection or thermal effects. The generated fractures improve the injectivity of the wellbore; however, they may also induce early water breakthrough at the production wells. If a fracture is created during water injection, it may remain open due to the decrease of the effective stress inside the fracture. When injection stops, the effective stress inside the fracture increases, causing fracture closure. Additionally, the injection of untreated water may cause pore-clogging at the surroundings of the injection point as minerals carried by the injection water precipitate. These combined effects (i.e., fracture closure and fracture face skin) can be captured using the proposed model.

Changqing field is a low-permeability reservoir in China. Since 1991 the field has been subjected to water injection to maintain reservoir pressure. Different studies carried out in the field indicate the presence of fracture creation and opening due to water injection (He et al., 2015). Fig. 20 (left) shows the pressure response of a fall-off test in a water injector well in the reservoir (Wang et al., 2018). The well injected water into a low permeability reservoir for some time. Later, the well was closed, and the pressure response was monitored. The measured pressure behavior points towards a well with water-induced fractures (WIF). The well and the adjacent wells were monitored, and it was reported

<figure><img src="image_24.png" /><figcaption>Fig. 19. Comparison of wellbore skin and fracture face skin on well performance.</figcaption></figure>

that there was no closed boundary near the well. The pressure and its derivative showed a storage behavior at early times, after which an upward in the pressure and its derivative was observed. The authors attribute this behavior to fracture conductivity reduction to fracture closure.

Using the proposed model, the pressure behavior can be matched for a fractured well going through fracture closure. Furthermore, the proposed model introduces the fracture face skin that water injector wells may experience due to the precipitation of the minerals or other formation damage. Fig. 20 (left) shows the matching of the field measurements and the solution obtained with the proposed model. The figures also show the results obtained by Wang et al. which model the effect of waterflood-induced fractures.

Table 4 lists the parameters obtained for pressure matching of the fall-off test. It is important to highlight that Wang et al. (2018), pointed out the high values of fracture-face skin ($S_f = 1$) obtained by a semi-analytical model for injector wells (Wang et al., 2017). Since the fracture-face skin usually ranges from 0 to 1 (Romero et al., 2003), in this case, a moderated value such as the one obtained with the model presented in this work ($S_f = 0.015$) could be more appropriate for a water injector well.

## 5.2. Horizontal hydraulic fractured well

Fig. 20 (right) shows the data and the matching for a horizontal well with four hydraulic fractures. The well was drilled in a shallow reservoir, which is a shelf-deposit of dolomite and anhydrite (Raghavan et al., 1997). The pressure derivative curve shows a short period of wellbore-storage effect, after which a lineal flow is identified by a slope of $1/2$, followed by a compound-linear flow at the late stage. Using the proposed model, the pressure behavior is matched, and Table 5 lists the parameters obtained for the pressure test

# 6. Conclusions

We developed a semi-analytical model to capture the pressure response from fractured wells exhibiting fracture closure and fracture face skin. The proposed model solves the diffusivity equation in the matrix formation using a semi-analytical approach based on superposition in space, allowing to include multistage fractures. A 1D numerical approach is used to solve the diffusivity equation within the hydraulic fractures. Pressure continuity is imposed through the fracture face skin (damage zone), allowing the crossflow between the reservoir matrix and the hydraulic fractures. Extensive testing with analytical and numerical simulations has been conducted to verify the proposed model.

Incorporating effective fracture compressibility and fracture conductivity reduction allows capturing the pressure response for fracture closure at different production times. Effective fracture compressibility has an important effect at early times in massive hydraulic fractures and compressible fractures, where the initial fracture flow cannot be considered incompressible. The same behavior was observed in high-storage fractures. The proposed model can reproduce the early response for a compressible fracture system, allowing its use for post-fracturing fracture characterization.

Unlike previous works, we showed that the fracture closure at high conductivity fractures is detectable when considering the minimum fracture conductivity $F_{CDmin}$, where the fracture reaches a finite conductivity behavior, and the pressure drop along the fracture cannot be neglected. For a low minimum fracture conductivity, the pressure response is evident even for infinite conductivity fractures, while at middle times, the pressure response can be better identified in a log-log plot. For high conductivity fractures or slow rates of closure, the pressure response at middle times is barely affected by fracture closure. We suggest using a semi-log approach to identify the rate and degree of fracture closure at late times. Furthermore, we proposed the use of an
<figure><img src="image_25.png" /><figcaption>**Fig. 20.** Type curve matching results for a water injector well with water-induced fractures (left), and type curve matching results for a multistage fractured well (right).</figcaption></figure>

<table><caption>**Table 4**
Type curve matching results for the field case.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Fracture face skin, S<sub>f</sub>[−]</td><td>0.015</td></tr><tr><td>Initial fracture conductivity, F<sub>CD</sub>[−]</td><td>80</td></tr><tr><td>Fracture conductivity at the end of the test, F<sub>CD</sub>[−]</td><td>4.2</td></tr><tr><td>Dimensionless permeability modulus, γ<sub>D</sub>[−]</td><td>3.2</td></tr><tr><td>Dimensionless fracture diffusivity, η<sub>D</sub>[−]</td><td>30,000</td></tr><tr><td>Dimensionless wellbore storage, C<sub>D</sub></td><td>7 × 10<sup>−4</sup></td></tr></tbody></table>

<table><caption>**Table 5**
Type curve matching results for the fractured horizontal well.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Fracture face skin in all fractures, S<sub>f</sub>[−]</td><td>0.01</td></tr><tr><td>Initial fracture conductivity, F<sub>CD</sub>[−]</td><td>15</td></tr><tr><td>Dimensionless permeability modulus, γ<sub>D</sub>[−]</td><td>2.8</td></tr><tr><td>Dimensionless fracture diffusivity, η<sub>D</sub>[−]</td><td>1 × 10<sup>5</sup></td></tr><tr><td>Dimensionless wellbore storage, C<sub>D</sub></td><td>0.03</td></tr></tbody></table>

effective wellbore radius once the radial flow is achieved. The pseudo-radial approach provides a fast and less computationally demanding method for characterizing the conductivity reduction at late times. The effective wellbore approach is only based on matching the pressure drop since the derivative in radial flow is constant. For the late times, the conductivity change is detectable for $F_{CDmin} < 50$.

The current models in the literature that use wellbore skin to mimic fracture face skin are inaccurate, and therefore, rigorous modeling is needed. Our model addresses this limitation and successfully incorporates the fracture face skin assuming a thin layer of reduced permeability surrounding the fracture walls. Additionally, the model can be used for horizontal multistage fracture wells, waterflood-induced fractured wells, and vertical fracture wells. Results showed that fracture face skin delays the start of the transition flow regime and introduces an additional pressure drop in the fracture. Consequently, an additional and earlier start of fracture conductivity reduction is triggered. The proposed PTA-based model offers a diagnostic tool to help identify and evaluate fracture closure and fracture face skin in hydraulically

fractured well.

The results from PTA test must be interpreted with caution based on the assumptions considered. The model assumes single-phase flow, which may not be valid in early flowback period after fracturing. The flowback flow is comprised by fracturing and reservoir fluid. Additionally, even though the reduction of permeability is accounted as an exponential decrease with respect to the effective stress, this should be confirmed by using laboratory tests or robust geomechanical models. Further, the presented model assumes a constant skin with time, which may not always be valid. A potential extension of this work is to consider the effect of multiphase flow and more complex consider boundary conditions.

### Credit author statement

**Ruben Figueroa Hernandez:** Writing – original draft preparation, Methodology, Visualization **Anis Younes:** Writing – review & editing, Validation **Marwan Fahs:** Methodology, Editing, Verification **Hussein Hoteit:** Conceptualization, Methodology, Writing – review & editing, Validation, Supervision.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgments

This publication is based upon work supported by the King Abdullah University of Science and Technology (KAUST) Research Funding Office under Award No. 4357.7he authors thank CMG Ltd. and KAPPA for offering the CMG STARS academic licenses, and Saphir software licenses, respectively, to KAUST.
# Appendix A. diffusivity equation in the matrix

The diffusivity equation in the matrix is given by:

$$ \frac{\partial^2 p_m}{\partial x^2} + \frac{\partial^2 p_m}{\partial y^2} = \frac{\varphi_m \mu c_{im}}{k_m} \frac{\partial p_m}{\partial t} \qquad (27) $$

To normalize the diffusivity equations, we introduce the following dimensionless groups, defined by.

## Dimensionless space and time variables

$$ x_D = \frac{x}{x_f}, y_D = \frac{y}{x_f}, t_D = \frac{k_m t}{\varphi_m \mu c_{im} x_f^2}, \qquad (28) $$

## Dimensionless pressures

$$ p_{mD} = \frac{2\pi k_m h(p_i - p_m)}{qB\mu}, \qquad (29) $$

The right-hand side of Eq. (27) can be expressed in terms of dimensionless variables. Using Leibniz rule, the partial derivatives can be expressed as:

$$ \frac{\partial p_m}{\partial t} = \frac{\partial p_m}{\partial p_D} \frac{\partial t_D}{\partial t} \frac{\partial p_{mD}}{\partial t_D} = - \frac{qB\mu}{2\pi k_m} \frac{k_m}{\mu \varphi_m c_{im} x_f^2} \frac{\partial p_{mD}}{\partial t_D}, \\ \frac{\partial p_m}{\partial x} = \frac{\partial p_f}{\partial x_D} \frac{\partial x_D}{\partial x} \frac{\partial p_{mD}}{\partial x_D} = - \frac{qB\mu}{2\pi k_m h} \frac{\partial p_{mD}}{\partial x_D}, \\ \frac{\partial^2 p_f}{\partial x^2} = - \frac{qB\mu}{2\pi k_m h} \frac{\partial^2 p_{mD}}{\partial x_D^2}. \qquad (30) $$

Using the definition of the dimensionless variables defined previously, Eq. (27) becomes,

$$ \frac{qB\mu}{2\pi k_m h x_f^2} \frac{\partial^2 p_{mD}}{\partial x_D^2} + \frac{qB\mu}{2\pi k_m h x_f^2} \frac{\partial^2 p_{mD}}{\partial y_D^2} = \frac{\varphi_m \mu c_{im}}{k_m} \frac{qB\mu}{2\pi k_m h x_f^2} \frac{k_m}{\varphi_m \mu c_{im}} \frac{\partial p_{mD}}{\partial t_D}, \qquad (31) $$

which simplifies to,

$$ \frac{\partial^2 p_{mD}}{\partial x_D^2} + \frac{\partial^2 p_{mD}}{\partial y_D^2} = \frac{\partial p_{mD}}{\partial t_D} \qquad (32) $$

# Appendix B. diffusivity equation in the HF

The fracture is assumed to be directly connected to the wellbore. We assume that the reservoir contributes flux only through the fracture, and there is no direct flow from the reservoir to the well. The fracture is modeled as a planar vertical fracture inside the reservoir. The fracture system contains geometric and hydraulic properties different from the matrix. To model the pressure behavior inside the fracture zone, we subdivide the fracture into control volumes, as displayed in Fig. A1.

<figure><img src="image_26.png" /><figcaption>Fig. A.1. Fracture system discretized into 1D control volumes within the HF in the x-direction, and connected to the skin zone in the y-direction.</figcaption></figure>

The control volume has a pore-volume of $PV = \varphi_f \Delta x \Delta y \Delta z$ in 3D. The fluid flows along the x-direction, and the matrix loads fluid to the HF through the skin zone in the y-direction (see Fig. A1). The fluid mass flux in the fracture ($q_f x$) across the interface of the control volume located at x is written as the fluid density ($\rho_x$) at x, multiplied by the interface area ($A_x = \Delta y \Delta z$) and the flow velocity ($v_f = - \frac{k_f}{\mu} \nabla p_f$) in the fracture, that is $q_f x = A_x \rho_x v_f x$. The mass conservation equation in the fracture's control volume in 3D can therefore be written as:

$$ (q_{f,x} - q_{f,x+\Delta x}) + (q_{s,y} - q_{s,y+\Delta y}) + (q_{f,z} - q_{f,z+\Delta z}) = \frac{\partial(\rho_f \Delta x \Delta y \Delta z)}{\partial t} \qquad (33) $$

With no flow occurring in the vertical z-direction, one gets ($q_{f,z} = q_{f,z+\Delta z} = 0$), and with the line-symmetry across the two sides of the HF, one gets ($q_{s,y+\Delta y} = - q_{s,y}$), which represents the matrix-to-fracture flux through the skin zone. Notice that the total flux $2q_s$ acts as a line-source (or sink) term to
the fracture. Then, Eq. (33) can be written as:

$$-(q_{x+\Delta x} - q_{f,x}) - 2q_{s,y+\Delta y} = \frac{\partial(\rho\varphi_f \Delta x \Delta y \Delta z)}{\partial t} \qquad (34)$$

Replacing the fluxes by their expressions in terms of velocity, one gets,

$$-(A_{x+\Delta x}\rho_{x+\Delta x}v_{f,x+\Delta x} - A_s\rho_s v_{f,x}) - 2A_{y+\Delta y}\rho_{y+\Delta y}v_{s,y+\Delta y} = \frac{\partial(\rho\varphi_f \Delta x \Delta y \Delta z)}{\partial t} \qquad (35)$$

In the above equation, $\Delta y$ represents the fracture aperture $w(\rho)$, which is a function of pressure (and therefore time), while $\Delta x$ and $\Delta z$ are independent of time and pressure. Therefore, dividing by $\Delta x \Delta y \Delta z$, Eq. (35) becomes in its continuous spatial form as,

$$-\frac{\partial(\rho_s v_{f,x})}{\partial x} - 2\frac{\rho_s v_{s,y}}{w} = \frac{1}{w} \frac{\partial(w\varphi_f \rho)}{\partial t} \qquad (36)$$

Using Darcy's velocity given by $v_{f,x} = -\frac{k_f}{\mu} \frac{\partial p_f}{\partial x}$, and $v_{s,y} = -\frac{k_s}{\mu} \frac{\partial p_s}{\partial y}$, Eq. (36) reads as,

$$\underbrace{\frac{\partial}{\partial x} \left( \rho \frac{k_f}{\mu} \frac{\partial p_f}{\partial x} \right) + \frac{2\rho}{w} \underbrace{\frac{k_s}{\mu} \frac{\partial p_s}{\partial y}}_{\text{matrix-to-fracture sink/source term}}}_{\text{sink/source term}} = \frac{1}{w} \frac{\partial(w\varphi_f \rho)}{\partial t} \qquad (37)$$

The accumulation term can be derived using Chain's rule as,

$$\begin{align}
\frac{\partial(w\varphi_f \rho)}{\partial t} &= \varphi_f \rho \frac{\partial w}{\partial p_f} \frac{\partial p_f}{\partial t} + \varphi_f w \frac{\partial \rho}{\partial p_f} \frac{\partial p_f}{\partial t} + w \rho \frac{\partial \varphi_f}{\partial p_f} \frac{\partial p_f}{\partial t}; \nonumber \\
&= \varphi_f \rho w \left( \frac{1}{w} \frac{\partial w}{\partial p_f} \right) \frac{\partial p_f}{\partial t} + \varphi_f \rho w \left( \frac{1}{\rho} \frac{\partial \rho}{\partial p_f} \right) \frac{\partial p_f}{\partial t} + \varphi_f \rho w \left( \frac{1}{\varphi_f} \frac{\partial \varphi_f}{\partial p_f} \right) \frac{\partial p_f}{\partial t}, \nonumber \\
&= \varphi_f \rho w \left( c_{\text{aperture}} + c_{\text{fluid}} + c_{\text{rock}} \right) \frac{\partial p_f}{\partial t}, \nonumber \\
&= \varphi_f \rho w c_f \frac{\partial p_f}{\partial t}. \tag{38}
\end{align}$$

Notice in the above equation, the total fracture compressibility is $c_{if} = (c_{aperture} + c_{fluid} + c_{rock})$, where $c_{rock}$, $c_{fluid}$, and $c_{aperture}$ are, respectively, rock formation compressibility in the fracture, fluid compressibility, and fracture compressibility reflecting the change in aperture as the pressure changes. As shown in the literature, fracture compressibility has a visible effect on the early stages of production (Liao et al., 2019) and a neglective effect in middle and late times. Assuming the early stage is relatively short compared with the total flow period, we assumed that the fracture compressibility is not a function of pressure. Additionally, to simplify the problem, we assume the fracture diffusivity is mainly governed by fracture compressibility, the permeability variation in the fracture diffusivity is not considered. Eq. (37) can then be linearized using the square-pressure, or pseudo-pressure approaches (see, e.g., (Lee, John, John B. Rollins, 2003)), such that:

$$\frac{\partial^2 p_f}{\partial x^2} + \frac{2}{w} \frac{k_s}{k_f} \frac{\partial p_s}{\partial y} = \frac{\varphi_f \mu c_{if}}{k_f} \frac{\partial p_f}{\partial t}. \qquad (39)$$

Because of the high permeability of the HF relative to the matrix, the pressure establishes the steady-state quickly within the fracture, which also reflects insignificant fluid accumulation within the fracture. Therefore, Eq. (39) reduces to,

$$\frac{\partial^2 p_f}{\partial x^2} = - \frac{2}{w} \frac{k_s}{k_f} \frac{\partial p_s}{\partial y}. \qquad (40)$$

To normalize the diffusivity equations, we introduce the following dimensionless groups, defined by:

Dimensionless pressures

$$p_{fD} = \frac{2\pi k h (p_i - p_f)}{q B \mu}, \qquad (41)$$

Dimensionless HF conductivity

$$F_{CD} = \frac{k_f w}{k_m x_f}.$$

We use Leibniz rule as in Eq. (30), the partial derivatives can be expressed as:

$$\frac{\partial p_s}{\partial y} = \frac{\partial p_s}{\partial y_D} \frac{\partial y_D}{\partial y} \frac{\partial p_{sD}}{\partial y_D} = \frac{q B \mu}{2 \pi k_h x_f} \frac{\partial p_{sD}}{\partial y_D} \qquad (42)$$

Substituting Eq. (42) in Eq. (39) and with the other terms as in Eq. (30), one gets,
$$ \frac{\partial^2 p_{jD}}{\partial x_D^2} + \frac{2}{F_{CD}} \frac{\partial p_{sD}}{\partial y_D} = \frac{1}{\eta_D} \frac{\partial p_{jD}}{\partial t_D}, \qquad (43) $$

where,

$$ \eta_D = \frac{k_f \varphi_{im} c_{im}}{k_m \varphi_f c_{if}}, \qquad (44) $$

## Appendix C. fracture face skin

Fracture face skin refers to a change in the reservoir permeability surrounding the fracture. This change in permeability is generated during proppant placement or scale precipitation. The damaged zone around the vertical fracture (see Fig. A1) is extended over a distance $b_s$ from the fracture face to the reservoir. This zone has a reduced permeability $k_s$, which generates an additional pressure drop $\Delta p_s$. In this model, the fracture face skin is accounted as an internal boundary condition for the matrix system.

<figure><img src="image_27.png" /><figcaption>Fig. C.1. Pressure profile versus distance from the HF showing the contribution of the skin zone to the pressure drop at the HF.</figcaption></figure>

To model the additional pressure drop within the skin zone, we assume that the damage zone is too narrow compared to the matrix size. Therefore, a steady-state pressure drop is assumed within the skin zone, i.e., between the matrix region and the HF region (see Fig. C1). The pressure drop caused by the damaged zone can be computed as:

$$ \Delta p_s = \Delta p_{real} - \Delta p_{ideal}, \qquad (45) $$

Where $\Delta p_{real}$, and $\Delta p_{ideal}$ are the actual and the skin-free pressure drops (see Fig. A2). Using Darcy's law to approximate this pressure drop, Eq. (45) becomes,

$$ \Delta p_s = \frac{qB\mu b_s}{4x_j h k_s} - \frac{qB\mu b_s}{4x_j h k_m} = \frac{qB\mu b_s}{4x_j h} \left( \frac{1}{k_s} - \frac{1}{k_m} \right) \qquad (46) $$

Note that in the above equation, $\Delta p_s$ is negative, zero, and positive when the permeability in the skin zone is larger, equal, and smaller than the matrix permeability, respectively. Consistently with the dimensionless pressure definition, we define the fracture face skin as.

$$ S_f = \frac{2\pi k_m h}{qB\mu} \Delta p_s = \frac{2\pi k_m h}{qB\mu} \left( \frac{qB\mu b_s}{4x_j h} \left( \frac{1}{k_s} - \frac{1}{k_m} \right) \right) = \frac{\pi b_s}{2x_j} \left( \frac{k_m}{k_s} - 1 \right) \qquad (47) $$

The matrix governing equation is coupled with the fracture equation using the internal boundary conditions, and then a system of equations is generated to solve the pressure response in each fracture slab (including the wellbore) and the rate coming from the matrix to the fracture.

## Appendix D. solution method

Using a fracture discretization since the fracture is divided into segments, the finite difference form of Eq. (17) for each fracture segment is given by:

$$ T_{Di-1,i} \bar{p}_{jDi-1} - (T_{Di,i-1} + T_{Di,i+1} + \chi_{Di}) \bar{p}_{jDi} + T_{Di+1,i} \bar{p}_{jDi+1} + \omega_{Di} \bar{q}_{jDi} = 0, \qquad (48) $$

Where,

$$ T_{Di,j} = \frac{\lambda_{Di} \lambda_{Dj}}{\lambda_{Di} + \lambda_{Dj}}, \quad \lambda_{D\theta} = \frac{F_{CD\theta}}{\Delta L_{jD\theta}/2}, \qquad (49) $$

$$ \chi_{D\theta} = s_D \left( \frac{F_{CD\theta} \Delta L_{jD\theta}}{\eta_{D\theta}} \right), \qquad (50) $$
$$ \omega_{D\theta} = -\pi \Delta L_{fD\theta} \qquad (51) $$

Eq. (48) requires special handling to account for the tips of the fracture, fracture interconnections, and the fracture segments connected to the wellbore. The discretized equation for the segment at the tip of the fracture reduces to:

$$ T_{Di-1,i} \bar{p}_{fDi-1} - (T_{Di,i-1} + \chi_{Di}) \bar{p}_{fDi} + \omega_{Di} \bar{q}_{fDi} = 0 \qquad (52) $$

For the fracture segments connected to the wellbore, the dimensionless bottom hole pressure $\bar{p}_{wD}$ is introduced to Eq. (48), leading to,

$$ \lambda_{Di} \bar{p}_{wD} - (\lambda_{Di} + T_{Di,i+1} + \chi_{Di}) \bar{p}_{fDi} + T_{Di+1,i} \bar{p}_{fDi+1} + \omega_{Di} \bar{q}_{fDi} = 0 \qquad (53) $$

For interconnecting fractures, the star-delta formation or other methods can be used to connect multiple fractures (Hoteit and Firoozabadi, 2005; Jiang et al., 2021; Karimi-Fard et al., 2004).

The total production rate from the well corresponds to the flux sum of all fracture blocks intersecting the well. For a well intersected by multiple fractures, Eq. (18) is modified to:

$$ \sum_{i \in \text{wellbore}} \frac{\partial \bar{p}_{fD}}{\partial x_D}_{xP=0,i} = - \frac{\pi}{F_{CD}^* s_D} \qquad (54) $$

The pressure response is obtained by solving the matrix-fracture system of equations. The matrix form of the discretized Eq. (54) is given by:

$$ \mathbf{B} \bar{p}_{fD} + \mathbf{C} \bar{p}_{wD} = \mathbf{D}, \qquad (55) $$

where $\mathbf{B}$ is a vector of dimension $N \times M$ which the $i$th element is $\lambda_i$ for the fracture segments connected to the wellbore. $\mathbf{C}$ is a vector of dimension 1 with elements equals to $- \sum_{i \in \text{wellbore}} \lambda_i$, and $\mathbf{D}$ is a vector of dimension 1 with element equals to $2\pi/s_D$. The final linear system is given in Eq. (22), such that,

$$ \begin{bmatrix} \mathbf{A} & -\mathbf{I} & \bar{\mathbf{0}} \\ \mathbf{F} & \mathbf{E} & \mathbf{B} \\ \bar{\mathbf{0}} & \mathbf{B}^{\prime} & \mathbf{C} \end{bmatrix} \begin{bmatrix} \bar{\mathbf{q}}_{fD} \\ \bar{\mathbf{p}}_{fD} \\ \bar{\mathbf{p}}_{wD} \end{bmatrix} = \begin{bmatrix} \bar{\mathbf{0}} \\ \bar{\mathbf{0}} \\ \mathbf{D} \end{bmatrix} $$

Where $\mathbf{F}$ is a diagonal matrix with dimensions ($N \times 2M, N \times 2M$) and diagonal value corresponds to $\omega_{Di}$ value Fig. D1 shows an example of the structure of the linear system corresponding to two HF with 4 fracture blocks.

<figure><img src="image_28.png" /><figcaption>Fig. D.1. Structure of the nonlinear system corresponding to two HF with four fracture blocks, where “*” indicates nonzero elements in the system.</figcaption></figure>

As the terms in the resulting system of equations are pressure-dependent $f(\delta_D)$, the generated system of equations is solved using an iterative procedure (Jiang et al., 2019b), where $\delta_D$ is defined in Eq. (10). The problem is linearized using the previous $\delta_D$ value for each fracture segment. An initial fracture conductivity $F_{CD}^* = F_{CDi}$ corresponds to a $\delta_D = 1$, and this first value is used to initialize the problem. The solution to the system of equations is computed assuming a $\bar{\delta}^k$ value, then using the computed dimensionless pressure, a new $\bar{\delta}^{k+1}$ is computed and compared with the previous $\bar{\delta}^k$, where $\bar{\delta}$ is a vector containing $\delta_D$ for each fracture segment. The iteration process finishes when $\bar{\delta}^k$ converges ($|\bar{\delta}^k - \bar{\delta}^{k+1}| < tol$). Once the pressures and rates are obtained in the Laplace space, the Stehfest algorithm converts the solutions to time-space. This procedure has to be repeated for each time step. Fig. D2 summarizes the iterative procedure to solve the system of equations.

Jiang et al. (2019b) performed a numerical analysis on the iterative procedure resulting in a recommendation of a tolerance criterion of $10^{-5}$ which was used in this work. Additionally, they suggested more numerical considerations for the computation in the iterative procedure.
<figure><img src="image_29.png" /><figcaption>Fig. D.2. Flow diagram for the iterative solution of the nonlinear system of equations.</figcaption></figure>

## References

Adegbola, K., Boney, C., 2002. Effect of fracture face damage on well productivity. In: International Symposium and Exhibition on Formation Damage Control. OnePetro.

Aguilera, R., 1999. Recovery factors and reserves in naturally fractured reservoirs. J. Can. Pet. Technol. 38.

Ahmed, U., Abou-Sayed, A.S., Jones, A.H., 1979. Experimental evaluation of fracturing fluid interaction with tight reservoir rocks and propped fractures. In: Symposium on Low Permeability Gas Reservoirs. OnePetro.

Al-Hussainy, R., Ramey, H.J., Crawford, P.B., 1966. The flow of real gases through porous media. J. Petrol. Technol. 18, 624–636.

Albattat, R., Hoteit, H., 2020. A Semi-analytical Approach to Model Drilling Fluid Leakage into Fractured Formation.

Albattat, R., He, X., AlSinan, M., Kwak, H., Hoteit, H., 2022. Modeling lost-circulation in fractured media using physics-based machine learning. In: 83rd EAGE Annu. Conf. Exhib., 1. https://doi.org/10.3997/2214-4609.202210204.-5.

CMG, 2019. CMOST: Intelligent Optimization & Analysis Tool.

Behr, A., Mtchedlishvili, G., Friedel, T., Haefner, F., 2006. Consideration of damaged zone in a tight gas reservoir model with a hydraulically fractured well. SPE Prod. Oper. 21, 206–211.

Bello, R., Wattenbarger, R., 2010a. Modelling and Analysis of Shale Gas Production with a Skin Effect. Journal of Canadian Petroleum Technology, OnePetro, pp. 37–48. https://doi.org/10.2118/143229-PA.

Bello, R., Wattenbarger, R., 2010b. Multi-Stage Hydraulically Fractured Shale Gas Rate Transient Analysis. North Africa Tech. https://doi.org/10.2118/126754-MS. Conf. Exhib.

Bennett, L., Calvez, J.L., Sarver, D.R.R., Tanner, K., Birk, W.S., Waters, G., Drew, J., Michaud, G., Primiero, P., Eisner, L., 2005. The source for hydraulic fracture characterization. Oilfield Rev. 17, 42–57.

Berumen, S., Tiab, D., 1997. Interpretation of stress damage on fracture conductivity. J. Pet. Sci. Eng. 17, 71–85.

Cardona, A., Finkbeiner, T., Santamarina, J.C., 2021. Natural rock fractures: from aperture to fluid flow. Rock Mech. Rock Eng. 54, 5827–5844. https://doi.org/10.1007/S00603-021-02565-1/FIGURES/3.

Charlez, P., Lemonnier, P., Ruffet, C., Bouteca, M.J., Tan, C., 1996. Thermally induced fracturing: analysis of a field case in North Sea. In: European Petroleum Conference. OnePetro.

Chen, H.Y., Poston, S.W., Raghavan, R., 1991. An application of the product solution principle for instantaneous source and green's functions. SPE Form. Eval. 6, 161–167.

Chen, D., Pan, Z., Ye, Z., 2015a. Dependence of gas shale fracture permeability on effective stress and reservoir pressure: model match and insights. Fuel 139, 383–392.

Chen, Z., Liao, X., Zhao, X., Dou, X., Zhu, L., 2015b. Performance of horizontal wells with fracture networks in shale gas formation. J. Pet. Sci. Eng. 133, 646–664.

Chen, Z., Liao, X., Yu, W., Zhao, X., 2018. Transient flow analysis in flowback period for shale reservoirs with complex fracture networks. J. Pet. Sci. Eng. 170, 721–737.

Cinco, H., Samaniego, V.F., 1981. Transient pressure analysis: finite conductivity fracture case versus damaged fracture case. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, San Antonio, Texas, US, pp. 4–7.

Cinco-Ley, H., 1981a. Transient pressure analysis: finite conductivity fracture case versus damaged fracture case. In: SPE Annual Technical Conference and Exhibition. OnePetro.

Cinco-Ley, H., 1981b. Transient pressure analysis for fractured wells. J. Petrol. Technol. 33, 1749–1766.

Cinco-Ley, H., Meng, H.-Z., 1988. Pressure transient analysis of wells with finite conductivity vertical fractures in double porosity reservoirs. In: SPE Annual Technical Conference and Exhibition. OnePetro.

Cinco-Ley, H., Samaniego, V.F., 1977. Effect of wellbore storage and damage on the transient pressure behavior of vertically fractured wells. In: SPE Annual Fall Technical Conference and Exhibition. Society of Petroleum Engineers.

Cipolla, C.L., Lolon, E.P., Erdle, J.C., Rubin, B., 2010. Reservoir modeling in shale-gas reservoirs. SPE Reservoir Eval. Eng. 13, 638–653.

Clarkson, C.R., 2013. Production data analysis of unconventional gas wells: review of theory and best practices. Int. J. Coal Geol. 109, 101–146.

Clarkson, C.R., Nobakht, M., Kaviani, D., Ertekin, T., 2012. Production analysis of tight-gas and shale-gas reservoirs using the dynamic-slippage concept. SPE J. 17, 230–242.

Clarkson, C.R., Qanbari, F., Nobakht, M., Heffner, L., 2013. Incorporating geomechanical and dynamic hydraulic-fracture-property changes into rate-transient analysis: example from the haynesville shale. SPE Reservoir Eval. Eng. 16, 303–316.

Clavijo, S.P., Addassi, M., Finkbeiner, T., Hoteit, H., 2022. A coupled phase-field and reactive-transport framework for fracture propagation in poroelastic media. Sci. Rep. 12, 1–15.

Cramer, D.D., 2005. Fracture skin: a primary cause of stimulation ineffectiveness in gas wells. In: SPE Annu. Tech. Conf. Exhib., https://doi.org/10.2118/96869-MS.

Ezulike, D.O., Dehghanpour, H., 2014. A model for simultaneous matrix depletion into natural and hydraulic fracture networks. J. Nat. Gas Sci. Eng. 16, 57–69.

Fan, L., Martin, R., Thompson, J., Atwood, K., Robinson, J., Lindsay, G., 2011. An integrated approach for understanding oil and gas reserves potential in eagle ford shale formation. In: Soc. Pet. Eng. - Can. Unconv. Resour. Conf. 2011, CURC 2011, 2, pp. 970–984. https://doi.org/10.2118/148751-MS.
Fu, P., Settgast, R.R., Hao, Y., Morris, J.P., Ryerson, F.J., 2017. The influence of hydraulic fracturing on carbon storage performance. J. Geophys. Res. Solid Earth 122, 9931-9949.

Gdanski, R.D., Weaver, J.D., Slabaugh, B.F., Walters, H.G., Parker, M.A., 2005. Fracture Face Damage-It Matters. In: SPE European Formation Damage Conference. OnePetro.

Gringarten, A.C., Ramey, H.J., Raghavan, R., 1974. Unsteady-state pressure distributions created by a well with a single infinite-conductivity vertical fracture. Soc. Petrol. Eng. J. 14, 347-360.

Guo, L., Fahs, M., Hoteit, H., Gao, R., Shao, Q., 2021. Uncertainty analysis of seepage-induced consolidation in a fractured porous medium. C. - Comput. Model. Eng. Sci. 129, 279-297. https://doi.org/10.32604/CMES.2021.016619.

Gupta, J.K., Zielonka, M.G., Albert, R.A., El-Rabaa, W., Burnham, H.A., Choi, N.H., 2012. Integrated methodology for optimizing development of unconventional gas resources. In: SPE Hydraulic Fracturing Technology Conference. OnePetro.

Hajirezaie, S., Wu, X., Peters, C.A., 2017. Scale formation in porous media and its impact on reservoir performance during water flooding. J. Nat. Gas Sci. Eng. 39, 188-202.

Hao, S., Kirk, M., Adwait, C., Dengen, Z., Hussein, H., 2015. Integrated Haynesville production analysis. In: Soc. Pet. Eng. - Unconv. Resour. Technol. Conf. URTEC. https://doi.org/10.2118/178722-ms, 2015.

Hawkes, R.V., Anderson, I., Bachman, R.C., Settari, A., 2013. Interpretation of Closure Stress in the Montney Shale Using PTA Based Techniques. https://doi.org/10.2118/163825-MS.

He, Y., Cheng, S., Li, L., Mu, G., Zhang, T., Xu, H., Qin, J., Yu, H., 2015. Waterflow direction and front characterization with multiple methods: a case study in Changqing oilfield, China. In: SPE Oil & Gas India Conference and Exhibition. OnePetro.

He, Y., Qin, J., Cheng, S., Chen, J., 2020. Estimation of fracture production and water breakthrough locations of multistage fractured horizontal wells combining pressure-transient analysis and electrical resistance tomography. J. Pet. Sci. Eng. 194, 107479.

He, X., Hoteit, H., AlSinan, M., Kwak, H.T., 2020a. Modeling hydraulic response of rock fractures under effective normal stress. In: Paper Presented at the ARMA/DGS/SEG International Geomechanics Symposium. SPE.

He, X., Hoteit, H., AlSinan, M.M., Kwak, H.T., 2020b. Modeling hydraulic response of rock fractures under effective normal stress. In: ARMA/DGS/SEG International Geomechanics Symposium 2020. IGS 2020.

He, X., Santoso, R., Zhu, W., Hoteit, H., AlSinan, M.M., Kwak, H.T., 2021. Coupled flow-normal-shear influence on fracture permeability: analysis and modeling. In: ARMA/DGS/SEG 2nd Int. Geomech. Symp.

He, Xupeng, Santoso, R., AlSinan, M., Kwak, H., Hoteit, H., 2021a. Constructing dual POROSITY models from high-resolution discrete-fracture models using deep neural networks. In: Society of Petroleum Engineers - SPE Reservoir Simulation Conference 2021. https://doi.org/10.2118/203901-MS. RSC 2021. OnePetro.

He, Xupeng, Sinan, M., Kwak, H., Hoteit, H., 2021b. A corrected cubic law for single-phase laminar flow through rough-walled fractures. Adv. Water Resour. 154, 103984 https://doi.org/10.1016/j.ADVWATRES.2021.103984.

Hoteit, H., Firoozabadi, A., 2005. Multicomponent fluid flow by discontinuous Galerkin and mixed methods in unfractured and fractured media. Water Resour. Res. 41, 1-15. https://doi.org/10.1029/2005WR004339.

Houze, O.P., Horne, R.N., Ramey, H.J., 1988. Pressure-transient response of an infinite-conductivity vertical fracture in a reservoir with double POROSITY behavior. SPE Form. Eval. 3, 510-518.

Huang, Y., Cheng, S., Yu, H., He, Y., Lin, B., Feng, N., 2017. A semi-analytical approach to estimate fracture closure and formation damage of vertically fractured wells in tight gas reservoir. J. Pet. Sci. Eng. 150, 85-90.

Hwang, Y., Lin, J., Schechter, D., 2013. Predicting Fractured Well Performance in Naturally Fractured Reservoirs. In: IPTC 2013: International Petroleum Technology Conference. European Association of Geoscientists & Engineers, p. 350.

Jia, P., Cheng, L., Clarkson, C.R., Qanbari, F., Huang, S., Cao, R., 2017. A Laplace-domain hybrid model for representing flow behavior of multifractured horizontal wells communicating through secondary fractures in unconventional reservoirs. SPE J. 22, 1856-1876.

Jia, P., Cheng, L., Clarkson, C.R., Huang, S., Wu, Y., Williams-Kovacs, J.D., 2018. A novel method for interpreting water data during flowback and early-time production of multi-fractured horizontal wells in shale reservoirs. Int. J. Coal Geol. 200, 186-198.

Jiang, L., Liu, T., Yang, D., 2019a. Effect of stress-sensitive fracture conductivity on transient pressure behavior for a horizontal well with multistage fractures. SPE J. 24, 1342-1363.

Jiang, L., Liu, T., Yang, D., 2019b. A semi-analytical model for predicting transient pressure behavior of a hydraulically fractured horizontal well in a naturally fractured reservoir with non-Darcy flow and stress-sensitive permeability effects. SPE J. 24, 1322-1341.

Jiang, L., Liu, J., Liu, T., Yang, D., 2020. Semi-analytical modeling of transient pressure behaviour for a fractured vertical well with hydraulic/natural fracture networks by considering stress-sensitive effect. J. Nat. Gas Sci. Eng. 82, 103477.

Jiang, L., Liu, J., Liu, T., Yang, D., 2021. Semi-analytical modeling of transient pressure behaviour for fractured horizontal wells in a tight formation with fractal-like discrete fracture network. J. Pet. Sci. Eng. 197, 107937.

Jones, F.O., 1975. A laboratory study of the effects of confining pressure on fracture flow and storage capacity in carbonate rocks. J. Petrol. Technol. 27, 21-27.

Karimi-Fard, M., Durlofsky, L.J., Aziz, K., 2004. An efficient discrete-fracture model applicable for general-purpose reservoir simulators. SPE J. 9, 227-236. https://doi.org/10.2118/88812-PA.

King, G.E., 2012. Hydraulic fracturing 101: what every representative, environmentalist, regulator, reporter, investor, university researcher, neighbor, and engineer should

know about hydraulic fracturing risk. J. Petrol. Technol. 64, 34-42. https://doi.org/10.2118/0412-0034-JPT.

Lawal, H., Jackson, G., Abolo, N., Flores, C., 2013. A novel approach to modeling and forecasting frac hits in shale gas wells. In: EAGE Annual Conference & Exhibition Incorporating SPE Europec. OnePetro.

Lee, John, Rollins, John B., J.P.S., 2003. Pressure Transient Testing.

Li, S., Cai, B., Ding, Y., Lu, Y., He, C., 2016. Damage of fracture face skin in massive hydraulic fracture stimulation. Soc. Pet. Eng. - IADC/SPE Asia pacific drill. In: Technol. Conf. https://doi.org/10.2118/180561-MS.

Liao, K., Zhang, S., Ma, X., Zou, Y., 2019. Numerical investigation of fracture compressibility and uncertainty on water-loss and production performance in tight oil reservoirs. Energies 12, 1189.

Lin, J., Zhu, D., 2014. Modeling well performance for fractured horizontal gas wells. J. Nat. Gas Sci. Eng. 18, 180-193.

Luo, L., Cheng, S., Lee, J., 2020. Characterization of refracture orientation in poorly propped fractured wells by pressure transient analysis: model, pitfall, and application. J. Nat. Gas Sci. Eng. 79, 103332.

Luo, L., Cheng, S., Lee, J., 2021. Semi-analytical model for characterization of anomalous linear-flow behavior in partially-propped fractured wells. J. Pet. Sci. Eng. 200, 108380.

Lynn, J.D., Nasr-El-Din, H.A., 1998. Evaluation of formation damage due to frac stimulation of a Saudi Arabian clastic reservoir. J. Pet. Sci. Eng. 21, 179-201.

Male, F., Islam, A.W., Patzek, T.W., Ikonnikova, S., Browning, J., Marder, M.P., 2014. Analysis of gas production from hydraulically fractured wells in the haynesville shale using scaling methods. In: Soc. Pet. Eng. - SPE USA Unconv. Resour. Conf., pp. 388-396. https://doi.org/10.2118/168993-MS, 2014.

Martins, J.P., Murray, L.R., Clifford, P.J., McLellan, W.G., Hanna, M.F., Sharp, J.W., 1995. Produced-water reinjection and fracturing in prudhoe bay. SPE Reservoir Eng. 10, 176-182.

Meng, M., Chen, Z., Liao, X., Wang, J., Shi, L., 2020. A well-testing method for parameter evaluation of multiple fractured horizontal wells with non-uniform fractures in shale oil reservoirs. Adv. Geo-energy Res. 4, 187-198.

Oelkers, E.H., Arkadaksiy, S., Afifi, A.M., Hoteit, H., Fedorik, J., Delaunay, A., Torres, J., E., Ahmed, Z.T., 2021. The subsurface carbonation potential of basaltic rocks from the Jizan region of Southwest Saudi Arabia. Hal Environ. Sci. 1-37.

Oelkers, E.H., Arkadaksiy, S., Afifi, A.M., Hoteit, H., Richards, M., Fedorik, J., Delaunay, A., Torres, J.E., Ahmed, Z.T., Kunnummal, N., Gislason, S.R., 2022. The subsurface carbonation potential of basaltic rocks from the Jizan region of Southwest Saudi Arabia. Int. J. Greenh. Gas Control 120, 103772. https://doi.org/10.1016/j.ijggc.2022.103772.

Omar, A.A., Addassi, M., Hoteit, H., Vahrenkamp, V., 2021. A new enhanced gas recovery scheme using carbonated water and supercritical CO₂. SSRN Electron 1-11. https://doi.org/10.2139/ssrn.3822140.

Ostensen, R.W., 1986. The effect of stress-dependent permeability on gas production and well testing. SPE Form. Eval. 1, 227-235.

Ozkan, E., Raghavan, R., 1991a. New solutions for well-test-analysis problems: part 1-analytical considerations (includes associated papers 28666 and 29213). SPE Form. Eval. 6, 359-368.

Ozkan, E., Raghavan, R., 1991b. New solutions for well-test-analysis problems: Part 2 computational considerations and applications. SPE Form. Evaluation 6, 369-378.

Palisch, T.T., Duenckel, R.J., Bazan, L.W., Heidt, J.H., Turk, G.A., 2007. Determining realistic fracture conductivity and understanding its impact on well performance-theory and field examples. In: SPE Hydraulic Fracturing Technology Conference. OnePetro.

Patzek, T.W., Male, F., Marder, M., 2013. Gas production in the Barnett Shale obeys a simple scaling theory. Proc. Natl. Acad. Sci. U.S.A. 110, 19731-19736. https://doi.org/10.1073/pnas.1313380110.

Pedrosa, O.A., 1986. Pressure transient response in stress-sensitive formations. In: SPE California Regional Meeting. OnePetro.

Prats, M., Hazebroek, P., Strickler, W.R., 1962. Effect of vertical fractures on reservoir behavior-compressible-fluid case. Soc. Petrol. Eng. J. 2, 87-94.

Qiao, T., Hoteit, H., Fahs, M., 2021. Semi-analytical solution to assess CO₂ leakage in the subsurface through abandoned wells. Energies 14, 2452.

Qin, J., Cheng, S., Li, P., He, Y., Lu, X., Yu, H., 2019. Interference well-test model for vertical well with double-segment fracture in a multi-well system. J. Pet. Sci. Eng. 183, 106412.

Raghavan, R.S., Chen, C.-C., Agarwal, B., 1997. An analysis of horizontal wells intercepted by multiple fractures. SPE J. 2, 235-245.

Reinicke, A., Rybacki, E., Stanchits, S., Huenges, E., Dresen, G., 2010. Hydraulic fracturing stimulation techniques and formation damage mechanisms—implications from laboratory testing of tight sandstone-proppant systems. Geochemistry 70, 107-117.

Romero, D.J., Valko, P.P., Economides, M.J., 2003. Optimization of the productivity index and the fracture geometry of a stimulated well with fracture face and choke skins. SPE Prod. Facil. 18, 57-64.

Roodhart, L.P., Kulper, T.O.H., Davies, D.R., 1988. Proppant-pack and formation impairment during gas-well hydraulic fracturing. SPE Prod. Eng. 3, 438-444.

Snæbjørnsdóttir, S.Ó., Sigfússon, B., Marieni, C., Goldberg, D., Gislason, S.R., Oelkers, E. H., 2020. Carbon dioxide storage through mineral carbonation. Nat. Rev. Earth Environ. 1, 90-102. https://doi.org/10.1038/s43017-019-0011-8.

Soliman, M.Y., 1986. Design and analysis of a fracture with changing conductivity. J. Can. Pet. Technol. 25.

Stehfest, H., 1970. Algorithm 368: numerical inversion of Laplace transforms [D5]. Commun. ACM 13, 47-49.

Sugar, A., Serag, M.F., Torrealba, V.A., Buttner, U., Habuchi, S., Hoteit, H., 2020. Visualization of polymer retention mechanisms in porous media UsingMicrofluidics.
In: Society of Petroleum Engineers - SPE Europec Featured at 82nd EAGE Conference and Exhibition. Society of Petroleum Engineers (SPE). https://doi.org/10.2118/200557-ms.

Sugar, A., Torrealba, V., Buttner, U., Hoteit, H., 2021. Assessment of polymer-induced clogging using microfluidics. SPE J. 26, 3793-3804. https://doi.org/10.2118/201626-PA.

Sun, H., Chawathé, A., Hoteit, H., Shi, X., Li, L., 2015. Understanding shale gas flow behavior using numerical simulation. SPE J. 20, 142-154. https://doi.org/10.2118/167753-PA.

Svendson, A.P., Wright, M.S., Clifford, P.J., Berry, P.J., 1991. Thermally induced fracturing of Ula water injectors. SPE Prod. Eng. 6, 384-390.

Tiab, D., Azzougen, A., Escobar, F.H., Bergumen, S., 1999. Analysis of pressure derivative data of finite-conductivity fractures by the Direct Synthesis technique. In: SPE Mid-continent Operations Symposium. OnePetro.

Tkhostov, B.A., Veziroba, A.D., Vendel'shteyn, B.Y., Dobrynin, V.M., 1970. Oil in Fractured Reservoirs, vols. 173-197. Izd. Nedra, Leningr.

Valdes-Perez, A.R., Pulido, H., Cinco-Ley, H., Larsen, L., 2013. A new double porosity fractal model for well test analysis with transient interporosity transference for petroleum and geothermal systems. In: Proceedings 38th Workshop on Geothermal Reservoir Engineering. Stanford University, Stanford, CA.

Van Everdingen, A.F., Hurst, W., 1949. The application of the Laplace transformation to flow problems in reservoirs. J. Petrol. Technol. 1, 305-324.

Verdon, J.P., Stork, A.L., 2016. Carbon capture and storage, geomechanics and induced seismic activity. J. Rock Mech. Geotech. Eng. 8, 928-935.

Wang, M., Leung, J.Y., 2015. Numerical investigation of fluid-loss mechanisms during hydraulic fracturing flow-back operations in tight reservoirs. J. Pet. Sci. Eng. 133, 85-102.

Wang, Y., Cheng, S., Feng, N., Xu, J., Qin, J., He, Y., Yu, H., 2017. Semi-analytical modeling for water injection well in tight reservoir considering the variation of waterflood-Induced fracture properties-Case studies in Changqing Oilfield, China. J. Pet. Sci. Eng. 159, 740-753.

Wang, Y., Cheng, S., Zhang, K., Xu, J., Qin, J., He, Y., Luo, L., Yu, H., 2018. Case studies: pressure-transient analysis for water injector with the influence of waterflood-induced fractures in tight reservoir. In: SPE Improved Oil Recovery Conference. OnePetro.

Wang, Y., Cheng, S., Zhang, K., Ayala, L.F., 2019. Investigation on the transient pressure response of water injector coupling the dynamic flow behaviors in the wellbore, waterflood-induced fracture and reservoir: semi-analytical modeling and a field case. Int. J. Heat Mass Tran. 130, 668-679. https://doi.org/10.1016/J.IJHEATMASTRANSFER.2018.09.083.

Wong, D.W., Harrington, A.G., Cinco-Ley, H., 1986. Application of the pressure derivative function in the pressure transient testing of fractured wells. SPE Form. Eval. 1, 470-480.

Wu, Z., Cui, C., Wang, Z., Sui, Y., Jia, P., Tang, W., 2018. Well testing model of multiple fractured horizontal well with consideration of stress-sensitivity and variable conductivity in tight gas reservoirs. Math. Probl Eng. 2018.

Xue, X., Yang, C., Sharma, V.K., Datta-Gupta, A., King, M.J., 2016. Reservoir and fracture flow characterization using a novel W(τ) formulation. SPE/AAPG/SEG unconv. In: Resour. Technol. Conf.. https://doi.org/10.15530/URTEC-2016-2440083, 2016.

Yang, D., Zhang, F., Styles, J.A., Gao, J., 2015. Performance evaluation of a horizontal well with multiple fractures by use of a slab-source function. SPE J. 20, 652-662.

Yu, W., Wu, K., Sepehrnoori, K., 2016. A semi-analytical model for production simulation from nonplanar hydraulic-fracture geometry in tight oil reservoirs. SPE J. 21, 1028-1040.

Zhang, F., Yang, D., 2018. Effects of non-Darcy flow and penetrating ratio on performance of horizontal wells with multiple fractures in a tight formation. J. Energy Resour. Technol. 140.

Zhang, Y., Yang, D., 2021. Evaluation of transient pressure responses of a hydraulically fractured horizontal well in a tight reservoir with an arbitrary shape by considering stress-sensitive effect. J. Pet. Sci. Eng. 202, 108518 https://doi.org/10.1016/J.PETROL.2021.108518.

Zhang, Z., He, S., Liu, G., Guo, X., Mo, S., 2014. Pressure buildup behavior of vertically fractured wells with stress-sensitive conductivity. J. Pet. Sci. Eng. 122, 48-55.

Zhou, W., Banerjee, R., Poe, B., Spath, J., Thambynayagam, M., 2014. Semi-analytical production simulation of complex hydraulic-fracture networks. SPE J. 19, 6-18.

Zoback, M.D., Byerlee, J.D., 1975. Permeability and effective stress. Am. Assoc. Petrol. Geol. Bull. 59, 154-158.