# Revising Transient-Pressure Solution for Vertical Well Intersected by a Partially Penetrating Fracture with Non-Darcy Flow Effect

Authors: Shuheng Cui; Jie Kong; Hongwei Yu; Cheng Chen; Junlei Wang
Journal: Geofluids
 (2020)

---

## Abstract

The principle purpose of this work is to formulate an accurate mathematical model to evaluate the transient pressure behavior of a well intercepted by a partially penetrating vertical fracture (PPVF) with non-Darcy flow effect. Fracture conductivity is taken into account by coupling the three-dimensional flow in reservoir and the two-dimensional flow within fracture; the Barree-Conway model is incorporated into the model to analyze non-Darcy flow behavior in fracture, which leads to the nonlinearity of the governing equations. A high-effective iterative algorithm using a combined technique of fracture-panel discretization and dimension transform is developed to render the nonlinear equations amenable to analytical linear treatment. On the basis of the solutions, the pressure response and its derivative type curves were generated to identify the evolution of flow regimes with time. Furthermore, the influences of fracture conductivity, penetration ratio, and non-Darcy characteristic parameters on pressure response are investigated. The results show that PPVF exhibits five typical flow regimes, and analytical solutions for each flow regime are similar to that for a fully penetrating vertical fracture (FPVF) that can be correlated with the penetration ratio and apparent conductivity. The non-Darcy flow effect is found to have more significant effect on the low and moderate conductivity, especially in early-stage flow regimes. When the penetration ratio is smaller than 0.5, the pressure behavior exhibit a more remarkable variation with penetration ratio. This study provides a better insight into understanding the influence of non-Darcy flow on flow regime identification.


---

## Body

## 1. Introduction

Fracturing stimulation has been widely applied in the development of tight hydrocarbon formations. Hydrocarbons are efficiently extracted from low-permeability formation into the wellbore throughout hydraulic fractures. It has been established that the inclination of the overall plane of a hydraulic fracture is determined by the axis of the least principal stress. If the axis is in the vertical or horizontal direction, the created fracture would be horizontal or vertical, respectively; otherwise, the fracture would be inclined [1]. In addition to those time-consuming numerical simulations, numerous works using various analytical or semianalytical methods have been published for study the flow behavior of vertical fracture [2–5], horizontal fracture [6–8], and inclined fracture [9–13]. These solutions serve as a theoretical basis of the pressure-transient and rate-transient analysis for hydraulically fractured well.At depths deeper than approximately 2000 ft in the unconventional reservoir, the direction of the least principal stress is generally parallel to the formation plane [14], and as a result, the hydraulic fractures are generated principally in the vertical direction. Although the intention was for oil companies to create fully penetrating hydraulic fractures, the actual measurement shows that hydraulic fracture may not extend throughout the entire vertical extent of the formation thickness [15, 16]. Only the effective height of the fracture that is propped open contributes to the production, and the three-dimensional flow pattern occurs in the reservoir. Alternatively, the flow in vertical direction should be considered. From the viewpoint of reservoir engineering, fracture dimensions (fracture length, width, and height) are all of much influence on the production performance. However, few studies have dealt with the effects of partially penetrating fracture. Raghavan et al. [17] presented a pioneering work of investigating the effect of fracture height of a single, vertical fracture. A uniform flux across the fracture plane was assumed, and an approximate analytical solution for the case of infinite conductivity was obtained by evaluating the uniform-flux solution at a special position. Rodriguez et al. [18] investigated the effect of the partial penetration of PPVF with infinite conductivity by discretizing the fracture into a set of uniform-flux planes and ignoring the pressure drop within fracture. Rodriguez et al. [19] further considered the pressure drop within finite conductivity fracture by establishing a Poisson’s type differential equation in fracture. Igbokoyi and Tiab [20] used the elliptical flow model in Laplace domain to obtain the solution for PPVF with infinite conductivity in naturally fractured reservoirs. Al-Rbeawi and Tiab successively used the type-curve matching technique and Tiab’s direct synthesis technique to analyze the pressure behavior of a horizontal well with multiple vertical partially hydraulic fractures [11, 12].Meanwhile, the non-Darcy flow effect should be considered because it usually happens due to the high-velocity flow within hydraulic fracture. As we know, non-Darcy flow within fracture has more important influence on transient pressure responses than non-Darcy flow in the reservoir [21, 22]. It is well demonstrated that non-Darcy effects account for a 35% decrease in productivity in a hydraulically fractured high-rate oil well, a larger productivity decrease for gas well [23]. It is well accepted that the fracture length and conductivity might be underestimated in the presence of non-Darcy flow. Umnuayponwiwat et al. [24] found that the errors of the estimation of fracture conductivity and fracture length would be up to 78% and 54% if neglecting non-Darcy flow effect for gas well. Therefore, it is of much practical importance to accurately evaluate the performance of appropriate conductivity fractured well with non-Darcy behavior. Forchheimer’s equation was widely used to analyze the flow behavior of finite conductivity FPVF with considering non-Darcy flow in many semianalytical approaches such as the work presented by Guppy et al. [25, 26]. Zeng and Zhao [27] used the Forchheimer number to quantify the effect of non-Darcy flow in the reservoir on the transient pressure behavior of vertical well through regarding it as rate-dependent skin factor. Valko and Amini [28] applied the 3D method of distributed volumetric sources to predict production from a horizontal well intersected by multiple transverse fractures with non-Darcy flow effect. Considered the fact that Forchheimer’s equation leads to an error at both low and high velocities [29, 30], and Zhang and Yang [31] incorporated the versatile Barree-Conway model to describe non-Darcy flow behavior in hydraulically fractured wells and developed an equation to correlate the Forchheimer equation and the Barree-Conway model. Jiang et al. [32] further applied it to quantify the non-Darcy flow effect in a hydraulically fracture horizontal well in a naturally fractured reservoir.These previous literatures demonstrated that there was a large reduction in the fracture conductivity when non-Darcy flow effect occurring in the FPVF was included. However, uncertainty still remains regarding whether the reduction becomes more significant as for the PPVF or not. It is necessary to emphasize that Zhang and Yang [33] used a novel slab source function to accurately quantify non-Darcy effect and penetrating ratio with consideration of the appropriate fracture dimension and its conductivity. The assumption of one-dimensional flow pattern in the finite conductivity fracture is proposed to simulate approximately the actual two-dimensional flow, which contributes to inaccurate results for the pressure-transient analysis of low/moderate-conductivity fracture.In this work, we develop an efficient and effective approach that provides accurate pressure-transient response of PPVF with non-Darcy flow effect during drawdown testing. Flux variations along the horizontal and vertical directions in the fracture are physically described by use of two-dimensional flow pattern. Different from the solutions using the Greens function method [11–13, 25, 31–33], a computational package based on Laplace solutions are derived to accelerate the calculation speed with high precision, and the technique using dimension transform is presented to eliminate the nonlinearity caused by non-Darcy effect. Subsequently, we identify the flow regimes and perform the sensitivity analysis of influence factors. In addition, a field case is used to illustrate the application of this model for accurate interpretation of well testing data.


## 2. Model Development

### 2.1. Model Illustrations

In this study, Figure1 shows the physical model of a vertical well intersected by a partially penetrating vertical fracture. There are several fundamental assumptions:
(i)
The formation is infinite in the lateral direction (x- and y-direction), and the boundaries on the top (z=h) and bottom (z=0) are impermeable. The formation is assumed to be homogeneous with constant porosity (φ), compressibility (ct), and permeability (km). The vertical permeability is set to be kmz(ii)
The fracture plane has a rectangular shape, and the compressibility of fracture is ignored(iii)
The length, width, height, and orientation angle of vertical fracture are denoted byLf, wf, hf, and θ, respectively. The starting coordinate is denoted by xofD,yofD,zofD.(iv)
The fluid is assumed to be single phase, isothermal, and slightly compressible(v)
The fluid flow is considered as Darcy’s law in the matrix, and the flow is described by non-Darcy flow in the fracture. Here, the non-Darcy flow is described with the Barree-Conway model [34](vi)
The well produces at a constant-rate conditionFigure 1
Schematic of a vertical well intersected by a partially penetrating fracture.Here, it is necessary to clarify the reason that we select the Barree-Conway model. First, the Forchheimer equation describes the laminar flow with inertia effect, which is given by(1)−∇p=μv1kfD⏟laminar+βρvμ⏟inertia=μv1kfapp,where β is non-Darcy flow coefficient, kfD represents the Darcy permeability, and kfapp represents the apparent permeability caused by non-Darcy flow effect. However, numerous studies demonstrate that 1/kfapp does not reduce to the inverse of the Darcy permeability 1/kfD when velocity approximates to zero, and kfapp approaches a value of zero at high velocities, rather than a limited value [29, 35]. This is inconsistent with Eq. (1), which indicates that Forchheimer’s equation may not be universal and only applicable in a limited range of velocities. Alternatively, the beta factor β is dependent on velocity, not a constant, and the kapp would deviate from the linear Forchheimer’s correlation at low and high velocities.The equation suggested by Barree and Conway [34] is introduced to recast the non-Darcy effect in the entire range of velocities. According to the Barree-Conway model, the ratio of apparent permeability to Darcy’s permeability is rewritten as follows:
(2)kfappkfD=kmr+1−kmr1+ρv/μτ=kmr+1−kmr1+FNDqcfD,where kmr is relative minimum permeability, τ is the characteristic length, and v is the fluid superficial velocity. In the form of dimensionless variables, the Reynolds number and the cross-sectional rate are defined, respectively, by
(3)FND=ρμτ×qrefwfh,qcfD=hhf×qcfqref.The Reynolds number is a constant, independent of fracture height. It is only determined by production rate and fracture width. Note that the fracture width is fixed in this study. Fracture conductivity is determined by changing the value of fracture permeability.


### 2.2. Mathematical Formulation

The fracture is represented by a rectangular porous medium of dimensionsLf×hf×wf. The flow parallel to the axis of fracture width can be ignored because of tiny fracture width, which signifies that the flux distribution can be assumed to be uniform in this direction [19]. After the non-Darcy effect is taken into account, according to the principle of mass balance, the fracture flow is described by the following partially differential equation in dimensionless form (dimensionless definitions are seen in Appendix A):
(4)∇kfappkfD∇pfD−2πCfDiqfDxfD,zfD+2πCfDiδxfD−xwfD∫0hfDqwfDzfD′δzfD−zfD′dzfD′=0.Laplace operator is written as∇=∂/∂xf+∂/∂zf.Initial dimensionless conductivity with regard to the Darcy permeability, which is a constant value, is defined by(5)CfDi=kfDwfkmLref.In Eq. (4), xwf is the wellbore location within fracture, qf is the flux-density function of spatial variables which represents the flux entering the fracture from the reservoir throughout the fracture surface. qwf is a source term representing the fluid extracted throughout the vertical wellbore within fracture, and it is parallel to the fracture height. In addition, the boundary conditions can be written as follows:
(6)∂pfD∂xfDxfD=xofD=∂pfD∂xfDxfD=xofD+LfD=0,∂pfD∂zfDzfD=zofD=∂pfD∂zfDzfD=zofD+hfD=0.In the vertical axis, the fracture plate is divided into a set of linear panels. Each panel is separated by red lines in Figure2. We eliminate the spatial dependence in zfD direction by integrating with respect to zfD varying from 0 to hfD:
(7)1hfD∫0hfD∇kappkD∇pfDdzfD−2πCfDihfD∫0hfDqfDxfD,zfDdzfD=−2πCfDihfDδxfD−xwfD∫0hfD∫0hfDqwfDzfD′δzfD−zfD′dzfD′dzfD.Figure 2
Discretization of PPVF along both the horizontal axis and the vertical axis.According to the studies provided by Al-Kobaisi et al. [36], the flow parallel to the fracture length dominates compared to the other direction. Therefore, the flux exchange between adjacent panels is ignored to some extent, i.e.,
(8)∂pfDn∂zfDzfD=zofDn=∂pfDn∂zfDzfD=zofDn+hfDn=0.As a result, thenth fracture panel forms the equation with regard to the weighted average of panel height hfDn, which is given by
(9)∂∂xfDnkappkD∂p¯fDn∂xfDn−2πCfDiq¯fDnxfDn+2πCfDinq¯wfDnδxfDn−xwfD=0.Equation (9) is an approximate treatment which converts the partially different equation for two-dimensional flow pattern into a system of ordinary-differential equations that is independent of zf-directional spatial variable. At a given position of xfD, the values of pressure in different panels might be different, which leads to a pressure difference contributing to the flux exchange in zf-direction. In fact, the approximation given in Eq. (9) has some inconsistencies, but it is close to the exact solution for the cases of interest; therefore, the equations could still be used for practical purposes. This approach is similar to the approximate analytical solution for the composite five-region model presented by Stalgorova and Mattar [37], and the following calculation results also verify our approximate approach.The solution for matrix system could be directly obtained by using the fundamental point-source solution in Laplace-transformed domain [38], which is given by
(10)p~DxD,yD,zD=∑n=1Nf∫0hfDn∫0LfDnq~DnξD,ζDK0srDndξDdζD+2∑n=1Nf∑m=1∞∫zofDnzofDn+hfDn∫0LfDnq~DnξD,ζDcosβmzDcosβmζDK0s+βm2hD2rDndξDdζD,where Nf is the number of fracture panels, βm=mπ, and the plane distance between the spatial point in the reservoir xD,yD and the spatial point on the fracture panel xofD+ξDcosθf,yofD+ξDsinθf is expressed as
(11)rDn=xD−xofD−ξDcosθf2+yD−yofD−ξDsinθf20.5.Although the Laplace-transformed solutions have great advantage in calculating the convolution, they usually pose computational problem.K0 is a zero-modified Bessel function of the second type. The case of rDn=0 leads to a singularity, and it is difficult to use numerical integration. Besides, since the integrals of Bessel functions appear in the terms of the infinite series [i.e., the second term on the right hand side of Eq. (10)], it is important to ensure the convergent and vanishing components of the infinite series. Some alterative solutions must be developed to accelerate the speed of the computations and improved the accuracy. An accurate computation package is provided in Appendix B, which is an important novelty of this work.


### 2.3. Semianalytical Solution

Although the two-dimensional equation is simplified into a system of one-dimensional equations, the governing equation for fracture panel described by Eq. (9) still has a strong nonlinearity nature because the term of kapp/kfD is a function of velocity with regard to temporal and spatial variables. We introduced the technique of dimension transformation presented by Luo and Tang [39] to render the nonlinear equation amenable to linear analytical treatment. Here, the definition of transformed dimension is given by
(12)ξDnxfDn=C^fDn,app⋅∫0xfDndxDCfDn,appxD,C^fDn,app=LfDn/∫0LfDndxDCfDn,appxD,and an apparent conductivity considering non-Darcy flow effect is defined as
(13)C^fDn,app=LfDn∫0LfDn1/CfDn,appxDdxD.Based on Eq. (2) and Eq. (5), the relation between apparent and initial dimensionless conductivity satisfies
(14)CfDn,appxfDn,qcfDn=CfDikmr+1−kmr1+FNDqcfDn.In this work, the conductivity based on the non-Darcy model is denoted asapparent conductivity, while the conductivity based on the Darcy model is denoted as (initial) conductivity. The apparent conductivity equals to the initial conductivity in the absence of non-Darcy flow (kmr=1orFND=0), and it will be smaller than the initial in the condition of non-Darcy flow (other cases).After substituting Eq. (12) into Eq. (9) and using Laplace transformation, a system of linear equations is obtained, which is
(15)∂2p¯~fDn∂ξDn2−2πC^fDn,appq¯~fDnξDn+2πC^fDn,appq¯~wfDnδξfDn−ξwfD=0.Eq. (15) is the well-known Fredholm integral equation. Integrating Eq. (15) with regard to ξDn from 0 to ξDn would yield the closed-form pressure solution and auxiliary solution (i.e., cross-sectional velocity), which are expressed, respectively, as follows:
(16)p¯~wD−p¯~fDnξDn=2πC^fDn,appq¯~wfDGξDn−ξwfDn−2πC^fDn,app∫ξwfDnξDndς∫0ςq¯~fDnζdζq¯~cfDnξDn=∫0ξDnq¯~fDnζdζ−q¯~wfDnHξDn−ξwfDn.Noting that according to the Darcy law in the transformed dimension, the dimensionless cross-sectional flow rate is defined asq¯~cfDn=hfDnC^fDn,app/2π∂p¯~fDn/∂ξDn.To obtain the unknowns, a semianalytical method is used to further discretize each panel intoNi segments with uniform flux (denoted by blue lines in Figure 2). According to the constant-rate condition, the constraint condition can be written as follows:
(17)∑n=1Nfq¯~wfDn=∑n=1Nf∑i=1Niq¯~fDn,iΔξDn,i=1s.In addition, the pressure and flow rate are continuous at the interface between the matrix and fracture. Here,(18)p¯~fDξDnxfDn,i=p~DxofD+xfDn,icosθ,yofD+xfDn,isinθ,zofDn+0.5hfDn.Because the solution of the fracture is derived from the new dimension, flux cannot be directly equalized at the interface. It is correlated as(19)q¯~fDn,i=q~Dn,i×ΔxfDn,iΔξDn,i.On the basis of the continuity condition and constraint condition, the matrix and fracture systems are coupled to generate closed-form equations. The unknown variables are listed as(i)
Ni×Nf uniform-flux rate of each segment, qfDn,i,i=1,2,⋯,Ni,n=1,2,⋯,Nf(ii)
Nf flow rate of each panel, qwfDn,n=1,2,⋯,Nf(iii)
The wellbore pressure,pwDThe total unknown variables areNi×Nf+Nf+1, which are presented in a vector form:
(20)XT=q¯~fD1,1,q¯~fD1,2,⋯,q¯~fD1,Ni,⋯,q¯~fDNf,Ni⏟Ni×Nf,q¯~wfD1,q¯~wfD2,⋯,q¯~wfDNf⏟Nf,p¯~wD.The structure of coefficient matrix and the known vector refers to our previous work [5].In each time step, an explicit iterative method is applied here to solve the nonlinear equations. Thus, Eq. (16) is written as
(21)p¯~wD−p¯~fDnκ+1ξDnκ=2πC^fDn,appκq¯~wfDκ+1GξDnκ−ξwfDnκ−2πC^fDn,appκ∫ξwfDnκξDnκdς∫0ςq¯~fDnκ+1ζdζq¯~cfDnκ+1ξDnκ=∫0ξDnκq¯~fDnκ+1ζdζ−q¯~wfDnκ+1HξDnκ−ξwfDnκ.The detailed process is illustrated as follows:(1)
Model inputs: reservoir parameters, fracture dimensions, and non-Darcy characteristic parameters(2)
Initialization: withκ=0, the apparent conductivity CfDn,appκ is assumed to be the initial conductivity CfDi. The q¯~fDnκ and q¯~cfDnκ are obtained by Gaussian elimination method and Stehfest numerical inversion(3)
Iterative process(a)
Calculating apparent conductivityCfDn,appκ and C^fDn,appκ and transforming xfDn into ξDnκ according to Eq. (12)(b)
Solving the explicit linear equations of Eq. (21) by Gaussian elimination method and Stehfest numerical inversion and obtaining the updated q¯~fDnκ+1, q¯~cfDnκ+1, and q¯~wfDκ+1(c)
Ifp¯~wDκ+1−p¯~wDκ<ε=10−5, then terminate the iterative process; otherwise, updating q¯~fDnκ=q¯~fDnκ+1 and q¯~cfDnκ=q¯~cfDnκ+1 with κ=κ+1 and return step (a) until convergence


## 2.1. Model Illustrations

In this study, Figure1 shows the physical model of a vertical well intersected by a partially penetrating vertical fracture. There are several fundamental assumptions:
(i)
The formation is infinite in the lateral direction (x- and y-direction), and the boundaries on the top (z=h) and bottom (z=0) are impermeable. The formation is assumed to be homogeneous with constant porosity (φ), compressibility (ct), and permeability (km). The vertical permeability is set to be kmz(ii)
The fracture plane has a rectangular shape, and the compressibility of fracture is ignored(iii)
The length, width, height, and orientation angle of vertical fracture are denoted byLf, wf, hf, and θ, respectively. The starting coordinate is denoted by xofD,yofD,zofD.(iv)
The fluid is assumed to be single phase, isothermal, and slightly compressible(v)
The fluid flow is considered as Darcy’s law in the matrix, and the flow is described by non-Darcy flow in the fracture. Here, the non-Darcy flow is described with the Barree-Conway model [34](vi)
The well produces at a constant-rate conditionFigure 1
Schematic of a vertical well intersected by a partially penetrating fracture.Here, it is necessary to clarify the reason that we select the Barree-Conway model. First, the Forchheimer equation describes the laminar flow with inertia effect, which is given by(1)−∇p=μv1kfD⏟laminar+βρvμ⏟inertia=μv1kfapp,where β is non-Darcy flow coefficient, kfD represents the Darcy permeability, and kfapp represents the apparent permeability caused by non-Darcy flow effect. However, numerous studies demonstrate that 1/kfapp does not reduce to the inverse of the Darcy permeability 1/kfD when velocity approximates to zero, and kfapp approaches a value of zero at high velocities, rather than a limited value [29, 35]. This is inconsistent with Eq. (1), which indicates that Forchheimer’s equation may not be universal and only applicable in a limited range of velocities. Alternatively, the beta factor β is dependent on velocity, not a constant, and the kapp would deviate from the linear Forchheimer’s correlation at low and high velocities.The equation suggested by Barree and Conway [34] is introduced to recast the non-Darcy effect in the entire range of velocities. According to the Barree-Conway model, the ratio of apparent permeability to Darcy’s permeability is rewritten as follows:
(2)kfappkfD=kmr+1−kmr1+ρv/μτ=kmr+1−kmr1+FNDqcfD,where kmr is relative minimum permeability, τ is the characteristic length, and v is the fluid superficial velocity. In the form of dimensionless variables, the Reynolds number and the cross-sectional rate are defined, respectively, by
(3)FND=ρμτ×qrefwfh,qcfD=hhf×qcfqref.The Reynolds number is a constant, independent of fracture height. It is only determined by production rate and fracture width. Note that the fracture width is fixed in this study. Fracture conductivity is determined by changing the value of fracture permeability.


## 2.2. Mathematical Formulation

The fracture is represented by a rectangular porous medium of dimensionsLf×hf×wf. The flow parallel to the axis of fracture width can be ignored because of tiny fracture width, which signifies that the flux distribution can be assumed to be uniform in this direction [19]. After the non-Darcy effect is taken into account, according to the principle of mass balance, the fracture flow is described by the following partially differential equation in dimensionless form (dimensionless definitions are seen in Appendix A):
(4)∇kfappkfD∇pfD−2πCfDiqfDxfD,zfD+2πCfDiδxfD−xwfD∫0hfDqwfDzfD′δzfD−zfD′dzfD′=0.Laplace operator is written as∇=∂/∂xf+∂/∂zf.Initial dimensionless conductivity with regard to the Darcy permeability, which is a constant value, is defined by(5)CfDi=kfDwfkmLref.In Eq. (4), xwf is the wellbore location within fracture, qf is the flux-density function of spatial variables which represents the flux entering the fracture from the reservoir throughout the fracture surface. qwf is a source term representing the fluid extracted throughout the vertical wellbore within fracture, and it is parallel to the fracture height. In addition, the boundary conditions can be written as follows:
(6)∂pfD∂xfDxfD=xofD=∂pfD∂xfDxfD=xofD+LfD=0,∂pfD∂zfDzfD=zofD=∂pfD∂zfDzfD=zofD+hfD=0.In the vertical axis, the fracture plate is divided into a set of linear panels. Each panel is separated by red lines in Figure2. We eliminate the spatial dependence in zfD direction by integrating with respect to zfD varying from 0 to hfD:
(7)1hfD∫0hfD∇kappkD∇pfDdzfD−2πCfDihfD∫0hfDqfDxfD,zfDdzfD=−2πCfDihfDδxfD−xwfD∫0hfD∫0hfDqwfDzfD′δzfD−zfD′dzfD′dzfD.Figure 2
Discretization of PPVF along both the horizontal axis and the vertical axis.According to the studies provided by Al-Kobaisi et al. [36], the flow parallel to the fracture length dominates compared to the other direction. Therefore, the flux exchange between adjacent panels is ignored to some extent, i.e.,
(8)∂pfDn∂zfDzfD=zofDn=∂pfDn∂zfDzfD=zofDn+hfDn=0.As a result, thenth fracture panel forms the equation with regard to the weighted average of panel height hfDn, which is given by
(9)∂∂xfDnkappkD∂p¯fDn∂xfDn−2πCfDiq¯fDnxfDn+2πCfDinq¯wfDnδxfDn−xwfD=0.Equation (9) is an approximate treatment which converts the partially different equation for two-dimensional flow pattern into a system of ordinary-differential equations that is independent of zf-directional spatial variable. At a given position of xfD, the values of pressure in different panels might be different, which leads to a pressure difference contributing to the flux exchange in zf-direction. In fact, the approximation given in Eq. (9) has some inconsistencies, but it is close to the exact solution for the cases of interest; therefore, the equations could still be used for practical purposes. This approach is similar to the approximate analytical solution for the composite five-region model presented by Stalgorova and Mattar [37], and the following calculation results also verify our approximate approach.The solution for matrix system could be directly obtained by using the fundamental point-source solution in Laplace-transformed domain [38], which is given by
(10)p~DxD,yD,zD=∑n=1Nf∫0hfDn∫0LfDnq~DnξD,ζDK0srDndξDdζD+2∑n=1Nf∑m=1∞∫zofDnzofDn+hfDn∫0LfDnq~DnξD,ζDcosβmzDcosβmζDK0s+βm2hD2rDndξDdζD,where Nf is the number of fracture panels, βm=mπ, and the plane distance between the spatial point in the reservoir xD,yD and the spatial point on the fracture panel xofD+ξDcosθf,yofD+ξDsinθf is expressed as
(11)rDn=xD−xofD−ξDcosθf2+yD−yofD−ξDsinθf20.5.Although the Laplace-transformed solutions have great advantage in calculating the convolution, they usually pose computational problem.K0 is a zero-modified Bessel function of the second type. The case of rDn=0 leads to a singularity, and it is difficult to use numerical integration. Besides, since the integrals of Bessel functions appear in the terms of the infinite series [i.e., the second term on the right hand side of Eq. (10)], it is important to ensure the convergent and vanishing components of the infinite series. Some alterative solutions must be developed to accelerate the speed of the computations and improved the accuracy. An accurate computation package is provided in Appendix B, which is an important novelty of this work.


## 2.3. Semianalytical Solution

Although the two-dimensional equation is simplified into a system of one-dimensional equations, the governing equation for fracture panel described by Eq. (9) still has a strong nonlinearity nature because the term of kapp/kfD is a function of velocity with regard to temporal and spatial variables. We introduced the technique of dimension transformation presented by Luo and Tang [39] to render the nonlinear equation amenable to linear analytical treatment. Here, the definition of transformed dimension is given by
(12)ξDnxfDn=C^fDn,app⋅∫0xfDndxDCfDn,appxD,C^fDn,app=LfDn/∫0LfDndxDCfDn,appxD,and an apparent conductivity considering non-Darcy flow effect is defined as
(13)C^fDn,app=LfDn∫0LfDn1/CfDn,appxDdxD.Based on Eq. (2) and Eq. (5), the relation between apparent and initial dimensionless conductivity satisfies
(14)CfDn,appxfDn,qcfDn=CfDikmr+1−kmr1+FNDqcfDn.In this work, the conductivity based on the non-Darcy model is denoted asapparent conductivity, while the conductivity based on the Darcy model is denoted as (initial) conductivity. The apparent conductivity equals to the initial conductivity in the absence of non-Darcy flow (kmr=1orFND=0), and it will be smaller than the initial in the condition of non-Darcy flow (other cases).After substituting Eq. (12) into Eq. (9) and using Laplace transformation, a system of linear equations is obtained, which is
(15)∂2p¯~fDn∂ξDn2−2πC^fDn,appq¯~fDnξDn+2πC^fDn,appq¯~wfDnδξfDn−ξwfD=0.Eq. (15) is the well-known Fredholm integral equation. Integrating Eq. (15) with regard to ξDn from 0 to ξDn would yield the closed-form pressure solution and auxiliary solution (i.e., cross-sectional velocity), which are expressed, respectively, as follows:
(16)p¯~wD−p¯~fDnξDn=2πC^fDn,appq¯~wfDGξDn−ξwfDn−2πC^fDn,app∫ξwfDnξDndς∫0ςq¯~fDnζdζq¯~cfDnξDn=∫0ξDnq¯~fDnζdζ−q¯~wfDnHξDn−ξwfDn.Noting that according to the Darcy law in the transformed dimension, the dimensionless cross-sectional flow rate is defined asq¯~cfDn=hfDnC^fDn,app/2π∂p¯~fDn/∂ξDn.To obtain the unknowns, a semianalytical method is used to further discretize each panel intoNi segments with uniform flux (denoted by blue lines in Figure 2). According to the constant-rate condition, the constraint condition can be written as follows:
(17)∑n=1Nfq¯~wfDn=∑n=1Nf∑i=1Niq¯~fDn,iΔξDn,i=1s.In addition, the pressure and flow rate are continuous at the interface between the matrix and fracture. Here,(18)p¯~fDξDnxfDn,i=p~DxofD+xfDn,icosθ,yofD+xfDn,isinθ,zofDn+0.5hfDn.Because the solution of the fracture is derived from the new dimension, flux cannot be directly equalized at the interface. It is correlated as(19)q¯~fDn,i=q~Dn,i×ΔxfDn,iΔξDn,i.On the basis of the continuity condition and constraint condition, the matrix and fracture systems are coupled to generate closed-form equations. The unknown variables are listed as(i)
Ni×Nf uniform-flux rate of each segment, qfDn,i,i=1,2,⋯,Ni,n=1,2,⋯,Nf(ii)
Nf flow rate of each panel, qwfDn,n=1,2,⋯,Nf(iii)
The wellbore pressure,pwDThe total unknown variables areNi×Nf+Nf+1, which are presented in a vector form:
(20)XT=q¯~fD1,1,q¯~fD1,2,⋯,q¯~fD1,Ni,⋯,q¯~fDNf,Ni⏟Ni×Nf,q¯~wfD1,q¯~wfD2,⋯,q¯~wfDNf⏟Nf,p¯~wD.The structure of coefficient matrix and the known vector refers to our previous work [5].In each time step, an explicit iterative method is applied here to solve the nonlinear equations. Thus, Eq. (16) is written as
(21)p¯~wD−p¯~fDnκ+1ξDnκ=2πC^fDn,appκq¯~wfDκ+1GξDnκ−ξwfDnκ−2πC^fDn,appκ∫ξwfDnκξDnκdς∫0ςq¯~fDnκ+1ζdζq¯~cfDnκ+1ξDnκ=∫0ξDnκq¯~fDnκ+1ζdζ−q¯~wfDnκ+1HξDnκ−ξwfDnκ.The detailed process is illustrated as follows:(1)
Model inputs: reservoir parameters, fracture dimensions, and non-Darcy characteristic parameters(2)
Initialization: withκ=0, the apparent conductivity CfDn,appκ is assumed to be the initial conductivity CfDi. The q¯~fDnκ and q¯~cfDnκ are obtained by Gaussian elimination method and Stehfest numerical inversion(3)
Iterative process(a)
Calculating apparent conductivityCfDn,appκ and C^fDn,appκ and transforming xfDn into ξDnκ according to Eq. (12)(b)
Solving the explicit linear equations of Eq. (21) by Gaussian elimination method and Stehfest numerical inversion and obtaining the updated q¯~fDnκ+1, q¯~cfDnκ+1, and q¯~wfDκ+1(c)
Ifp¯~wDκ+1−p¯~wDκ<ε=10−5, then terminate the iterative process; otherwise, updating q¯~fDnκ=q¯~fDnκ+1 and q¯~cfDnκ=q¯~cfDnκ+1 with κ=κ+1 and return step (a) until convergence


## 3. Results and Discussions

### 3.1. Model Validation

Before conducting the validation, a sensitivity analysis was first conducted to determine the number of fracture panels (Nf) and the number of fracture segments (Ni) needed to yield accurate results. As we know, both the number of fracture segment and the number of fracture panels increase with the decrease of conductivity. The following dimensionless fracture and reservoir data are used in this model: hD=2h/Lf=1, hfD=hf/h=1, km/kmz=1, CfDi=π, kmr=0.1, and FND=100, which indicate the case of low conductivity and strong non-Darcy effect. Figure 3 shows the dimensionless pressure from small to large time scope under different number of fracture panels. This indicates that the pressure drop approximate a stable and accurate result when the number of panels and segments is more than 1000 (10×100); therefore, the number of fracture panels of Nf=10 and the number of fracture segments Ni=100 are used in the following studies.Figure 3
Effect of the number of fracture panels on the calculation accuracy.In this section, existing solutions are used to verify our model. First, the results for finite-conductivity PPVF without non-Darcy flow effect were compared against results previously presented by Rodriguez et al. [19], where we set kmr=1 and FND=0 that assume the Darcy flow within fracture. Second, the results for finite-conductivity FPVF with non-Darcy flow effect using the Forchheimer equation were compared against the results presented by Guppy et al. [25] and Luo and Tang [39]; here, the Barree-Conway equation is simplified to the Forchheimer equation when kmr=0 and qDNDf=FND. Noting that fully penetrated fracture is the limiting case of our model. As seen in Figure 4, the calculated results from our model agree well with those published calculations for times of interest. The validations indicate that the proposed semianalytical model is reliable in simulating pressure transient response of finite-conductivity PPVF with non-Darcy effect.Figure 4
Comparison between the dimensionless pressures calculated by our solution and the alternative results proposed by (a) Guppy et al. [25] and (b) Luo and Tang [39].
(a)(b)In addition, the speed of the computations of two models is compared using the identical hardware platform. For the above simulations, the computation times using the improved algorithm in this paper and the previous algorithm are, respectively, 19.5 s and 28.8 s.


### 3.2. Transient-Flow Behavior

Figure5 shows the dimensionless pressure and derivative responses of the PPVF with a moderate conductivity (CfDi=100π) and a small penetration ratio (hfD=0.1) under a weak non-Darcy flow effect (kmr=0.5,FND=100). For the sake of comparison and description, the corresponding pressure responses caused by the Darcy flow are also presented. The black lines represent the non-Darcy case, while the red lines represent the Darcy case.Figure 5
Pressure and pressure derivative responses for PPVF with and without the non-Darcy flow effect.As analyzed in Figure5, the flow regimes are clearly distinguished by identifying the slopes on log-log curve of pressure-derivative in the Darcy condition. Five typical flow regimes are exhibited in sequence: (1) bilinear-flow period characterized by a 1/4-slope straight line, (2) formation-linear flow period identified by a 1/2-slope straight line, (3) early-radial flow regime, (4) compound-linear flow regime, and (5) pseudoradial flow regime. A detailed illustration was presented.


#### 3.2.1. Bilinear flow

The fluids beyond the fracture face in vertical direction are not developed as show in Figure6(a). An explicit solution for a FPVF was provided by Cinco-Ley and Samaniego [40] during bilinear flow:
(22)pwDtD=πΓ5/42CfDitD1/4.Figure 6
Schematic of some special flow regimes occurring during the production for a PPVF.
(a)
Bilinear and linear flow(b)
Early-radial flow(c)
Compound-linear flowFor PPVF with non-Darcy flow effect, the fracture behaves like FPVF with the same height (hf) during this period, so the solution is corrected as follows:
(23)hfDpwDtD=2πkmhfpi−pwqrefμ=πΓ5/42CfD,apptD1/4,where CfD,app is a correlation for the apparent fracture conductivity, which is the function of initial conductivity (CfDi), the Reynolds number (FND), and the minipermeability plateau (kmr). After simulating numerous cases, the relationship between apparent and initial (true) conductivities that is suitable for the special case using the parameters seen in Figure 5 is achieved, which is
(24)CfD,appCfDi=kmr+1−kmr1+2.51FND.


#### 3.2.2. Formation-linear flow

During this period as shown in Figure6(a), formation fluid enters the fracture in a direction perpendicular to the fracture surface, but the pressure response within fracture is not considered again. Note that the flow pattern in the fracture exhibits line-shaped during bilinear and linear flow regimes, which is confirmed in Figure 7(a). Based on the analytical solution for FPVF presented by Gil et al. [41], both penetration ratio and non-Darcy effect are further taken into account in this work. Similar to the derivation of bilinear-flow solution, we recast analytical solution for formation-linear-flow period as follows:
(25)hfDpwDtD=πtD+αCfD,app,where α is a constant given by the following: α=π/3 for CfD≥25, α=0.944 for 10≤CfD<25, and α=0.902 for 5≤CfD<10 [41].Figure 7
Pressure fields within two-dimensional fracture at (a)tD=10−6 and (b) tD=10−2.
(a)(b)


#### 3.2.3. Early-radial (ER) flow

If penetration ratio is small enough,hfD should be smaller than 0.5 suggested by Al-Rbeawi and Tiab [11]; an ER flow regime is developed around and along the fracture length in YZ plane. The ER flow regime is identified by a horizontal line on pressure derivative plot with
(26)tD×pwD′ER=0.5.During this period, fracture behaves like a finite-conductivity horizontal well with axis along the formation boundary, and the reservoir fluids flow radially inYZ plane towards fracture as seen in Figure 6(b). The effective penetration length becomes fracture length rather than formation thickness or fracture height, which is given by
(27)pwD,ER=pwD2hD=2πkmLfpi−pwqrefμ.The corresponding horizontal line on derivative curve is corrected as follows:(28)tD×pwD,ER′ER=0.52hD.Note that the flow pattern deviates from the line shaped, becoming the radial shaped as shown in Figure7(b).


#### 3.2.4. Compound-linear (CL) flow

CL flow often occurs after ER flow when both upper and bottom impermeable boundaries are felt. The fluid flow is predominantly perpendicular to the fracture face inXY plane, and the response does not reach the region beyond fracture tip. The PFVF behaves like a FFVF during this period, as shown in Figure 6(c). The characteristic of this flow behavior is very similar to the compound linear flow for multistage-fractured horizontal well presented by Chen and Raghavan [4] and fractured inclined well presented by Dinh and Tiab [10], which is identified by an approximated 1/2-slope straight line on log-log plot of pressure derivative.


#### 3.2.5. Pseudoradial (PR) flow

During this period, the reservoir flow begins to converge radially to the partial fracture. The flow inXY plane dominates compared with YZ plane; alternatively, the flow in z-direction can be ignored. The PR flow is characterized with tD×p’wD=0.5 on the log-log plot.Depending on the reservoir, fracture, and non-Darcy flow parameters used in a single test, the duration of some flow regime might elongate or shorten, even disappear. The non-Darcy flow effect leads to extra pressure drop, but cannot alter the characterization of flow regimes identified by the Darcy behavior.


### 3.3. Sensitivity Analysis

According to Eq. (2), the effect of non-Darcy flow is determined by two dimensionless parameters (i.e., kmr and FND) with the Barree-Conway model. It is assumed that the formation thickness, fracture length, and initial conductivity are given, i.e., hD=1, hfD=0.1, and CfDi=100π. To highlight the improvement of this work, the effect of dimensionless parameters with the Barree-Conway model and the penetration ratio of fracture are mainly discussed here because effects of dimensionless thickness, fracture conductivity, fracture location, and asymmetrical wellbore have been extensively investigated and analyzed [17, 25].Figure8 shows the influence of non-Darcy flow effect on transient pressure behavior of PFVF. Figure 8(a) displays the pressure and derivative curves for FND=0.1, 1, 10, and 100 at kmr=0.1. The Reynolds number (FND) is proportional to production rate according to Eq. (3). Increasing FND represents the increase of production rate. Besides, the case that FND=0 or kmr=1 indicates Darcy’s flow equation. The limiting Darcy cases are introduced, which are denoted by blue dashed lines. Case 5 indicates the pressure responses for the initial conductivity in the Darcy flow condition, while case 6 indicates the minimum conductivity in the Darcy flow condition (CfDmin=CfDi×kmr). For case 1~4, non-Darcy effect causes a significant reduction in fracture conductivity in simulated cases. As a result, an extra pressure drop is caused by non-Darcy effect. Case 4 exhibits a whole sequence of flow regimes similar to the initial case (case 5) due to a weak non-Darcy effect. The increasing of FND shortens the duration of formation-linear flow regime until disappears (i.e., case 1 and case 2). Meanwhile, the pressure drop and derivative curves also increase and approach the minimum case (case 6) because of strong non-Darcy effect. As a result, the curves of pressure responses with non-Darcy effect (i.e., case 1~case 4) are distributed within the range between the initial case (case 5) and the minimum case (case 6). Note that case 1 is almost overlapped with the minimum case for times of interest.Figure 8
Effect of the Reynolds number on (a) pressure response and (b) stabilized flux distribution.
(a)(b)The flux distribution stabilizes and remain unchanged when pseudoradial flow is established, which is shown in Figure8(b). Note that the flux distribution indicates the arithmetic average value of flux distribution among fracture panels. The integral with regard to xfD represents the volume of fluids entering the fracture. In case 6, most of the flow occurs from the half of fracture away from the wellbore. With the FND increasing, the characterization of fluid flow is approaching the unified case (case 6): the flow from the tips decreases, and more fluid comes from the half of fracture closer to the wellbore.Figure9(a) plots the curves of pressure responses for kmr=0.1, 0.2, 0.5, and 0.9 at FND=10. Different values of kmr represent different minimum cases under Darcy flow, which are denoted by blue dashed lines. The initial case is denoted by green dashed line. Solid lines indicate the cases with non-Darcy effect. The non-Darcy cases are located within the range between the initial case and the minimum case. As seen in Figure 9(a), as the value of kmr decreases, the pressure drop and its derivative would deviate from the initial case and approach the corresponding minimum case due to intense non-Darcy effect. When the value of kmr approaches unity, such as case 4, the pressure response of the minimum Darcy case is similar to the initial Darcy case. The apparent conductivity for case 4 is always in the order of magnitude for infinite conductivity during the process of conductivity degradation. As a result, the effect of non-Darcy flow tends to disappear. When the value of kmr is relatively smaller, the non-Darcy effect becomes significant. Figure 9(b) shows the stabilized average flux distribution along fracture. As the value of kmr decreases with a strong non-Darcy effect, the apparent conductivity decreases, and more fluid comes from the half of fracture away from the wellbore. This is different from the phenomenon described in Figure 8(b), because the minimum conductivity is different, not a unified value.Figure 9
Effect of the relative minimum permeability on (a) pressure response and (b) stabilized flux distribution.
(a)(b)As analyze from Figures8 and 9, the dimensionless parameters of larger FND and small kmr make the non-Darcy effect more significant. With the consideration of non-Darcy effect, the apparent fracture conductivity exhibits different conductivities in the Darcy flow condition when it actually has only one true value. Figure 10 further investigates the effect of penetration ratio in different conditions of initial conductivity. The value of FND and kmr is still the same as that in Figure 10, and the transient pressure and its derivative show a similar trend. For smaller conductivity (CfDi=10π) as shown in Figure 10(a), the pressure drop is always lower than the minimum Darcy case but higher than the initial Darcy case all the time in the condition of hfD=0.1. However, the non-Darcy case would overlap with the initial Darcy case in the late-time period in the condition of hfD=1. As a comparison, a higher conductivity is considered as shown in Figure 10(b). The pressure response of non-Darcy case deviates from the initial Darcy case, but would overlap with the initial case on the onset of ER flow regime in the condition of hfD=0.1. The pressure response of non-Darcy case might overlap with the initial case on advance in the condition of hfD=1.Figure 10
Effect of the penetration ratio on pressure response in the condition of (a) low and (b) high conductivity.
(a)(b)In summary, the effect of non-Darcy flow has a more significant influence on the condition of smaller penetration and lower initial conductivity. Alternatively, the effect of non-Darcy flow has a more significant influence on the condition of small penetration.


### 3.4. Application to Field Case

The case presents a pressure drawdown test in an oil well in the Tahe oil field of Xinjiang, China. The vertical well is hydraulically fractured, and the fracture is recognized as a partial penetrating vertical fracture according to the results of microseismic monitoring. The basic parameters used in this case are summarized in Table1.Table 1
Reservoir and fluid properties for the well.
ParameterValueUnitFormation thickness,h32.5mVolume factor,B1.05Formation porosity,φ19.8Fluid viscosity,μ2.0mPa•sWellbore radius,rw0.1mTotal compressibility,ct4.7×10−4MPa-1Production rate,qw10.0m3/dFigure11 provides the pressure drawdown and its derivative data on log-log plot. According to the flow regime identification in Section 3.2, it is shown from the derivative curve that a wellbore-storage effect is firstly recognized with a unity slope. As a result, the bilinear flow regime is masked by the wellbore-storage period. Then, the formation-linear flow regime appears, followed by a significant ER flow regime. The slope of ER flow is a negative value not zero, and the value of tD×p’wD is approximately 0.25. According to Eq. (28), the penetration ratio is about 1/4, so the fracture height is about 8.125 m height. Meanwhile, the negative slope for ER flow regime satisfies the characteristic of low and moderate conductivity of PPVF as presented in Figure 10(a). Eventually, the compound-linear flow and pseudoradial flow are found successively.Figure 11
Pressure transient analysis for (a) flow regime identification and (b) type-curve matching of the pressure drawdown data.
(a)(b)To match the entire drawdown data, the effects of wellbore storage capacity and fracture damage should be incorporated into our model. Because we use Laplace-domain solution, the solution accounting for these effects can be flexibly written as follows:(29)p~owD=sp~wD+Sds+CDs2sp~wD+Sd.Through constantly calculating the unknown variables, Figure11(b) presents the measured and final simulated pressure drawdown data, and an excellent agreement exists between the measured and simulated pressure for drawdown test. This case uses the responses of several different flow regimes at the same flow rate, instead of the responses of single flow regimes at two different rates as presented by Guppy et al. [25]. Nevertheless, the non-Darcy flow responses can be still matched with the Darcy flow response for a lower conductivity because of the nonuniqueness problem. If non-Darcy effect is not considered in the matching denoted by dashed lines, the dimensionless conductivity (CfDi) is calculated as 125.56.If non-Darcy effect is considered denoted by solid lines, the dimensionless initial conductivity (CfDi) is calculated as 204.29 with kmr=0.187, and FND=71.56 all the time. The resulting matching parameters are given in Table 2. Considered that the production rate of this well is relatively high, the fracture might be subject to the effects of non-Darcy flow. The traditional analysis method would result in nearly 40% error in the estimation of conductivity. The type curves based on the Darcy model would yield an underestimated conductivity.Table 2
Type-curve matching results for the drawdown test.
Without non-Darcy effectWith non-Darcy effectParameterValueUnitParameterValueUnitWellbore storage,C0.0012m3/MPaWellbore storage,C0.0024m3/MPaSkin factor,Sd0Skin factor,Sd0Fracture conductivity,kfwf65.99D•mFracture conductivity,kfwf284.86D•mFracture length,Lf58.4mFracture length,Lf66.4mFormation permeability,km0.018DFormation permeability,km0.042D


## 3.1. Model Validation

Before conducting the validation, a sensitivity analysis was first conducted to determine the number of fracture panels (Nf) and the number of fracture segments (Ni) needed to yield accurate results. As we know, both the number of fracture segment and the number of fracture panels increase with the decrease of conductivity. The following dimensionless fracture and reservoir data are used in this model: hD=2h/Lf=1, hfD=hf/h=1, km/kmz=1, CfDi=π, kmr=0.1, and FND=100, which indicate the case of low conductivity and strong non-Darcy effect. Figure 3 shows the dimensionless pressure from small to large time scope under different number of fracture panels. This indicates that the pressure drop approximate a stable and accurate result when the number of panels and segments is more than 1000 (10×100); therefore, the number of fracture panels of Nf=10 and the number of fracture segments Ni=100 are used in the following studies.Figure 3
Effect of the number of fracture panels on the calculation accuracy.In this section, existing solutions are used to verify our model. First, the results for finite-conductivity PPVF without non-Darcy flow effect were compared against results previously presented by Rodriguez et al. [19], where we set kmr=1 and FND=0 that assume the Darcy flow within fracture. Second, the results for finite-conductivity FPVF with non-Darcy flow effect using the Forchheimer equation were compared against the results presented by Guppy et al. [25] and Luo and Tang [39]; here, the Barree-Conway equation is simplified to the Forchheimer equation when kmr=0 and qDNDf=FND. Noting that fully penetrated fracture is the limiting case of our model. As seen in Figure 4, the calculated results from our model agree well with those published calculations for times of interest. The validations indicate that the proposed semianalytical model is reliable in simulating pressure transient response of finite-conductivity PPVF with non-Darcy effect.Figure 4
Comparison between the dimensionless pressures calculated by our solution and the alternative results proposed by (a) Guppy et al. [25] and (b) Luo and Tang [39].
(a)(b)In addition, the speed of the computations of two models is compared using the identical hardware platform. For the above simulations, the computation times using the improved algorithm in this paper and the previous algorithm are, respectively, 19.5 s and 28.8 s.


## 3.2. Transient-Flow Behavior

Figure5 shows the dimensionless pressure and derivative responses of the PPVF with a moderate conductivity (CfDi=100π) and a small penetration ratio (hfD=0.1) under a weak non-Darcy flow effect (kmr=0.5,FND=100). For the sake of comparison and description, the corresponding pressure responses caused by the Darcy flow are also presented. The black lines represent the non-Darcy case, while the red lines represent the Darcy case.Figure 5
Pressure and pressure derivative responses for PPVF with and without the non-Darcy flow effect.As analyzed in Figure5, the flow regimes are clearly distinguished by identifying the slopes on log-log curve of pressure-derivative in the Darcy condition. Five typical flow regimes are exhibited in sequence: (1) bilinear-flow period characterized by a 1/4-slope straight line, (2) formation-linear flow period identified by a 1/2-slope straight line, (3) early-radial flow regime, (4) compound-linear flow regime, and (5) pseudoradial flow regime. A detailed illustration was presented.


### 3.2.1. Bilinear flow

The fluids beyond the fracture face in vertical direction are not developed as show in Figure6(a). An explicit solution for a FPVF was provided by Cinco-Ley and Samaniego [40] during bilinear flow:
(22)pwDtD=πΓ5/42CfDitD1/4.Figure 6
Schematic of some special flow regimes occurring during the production for a PPVF.
(a)
Bilinear and linear flow(b)
Early-radial flow(c)
Compound-linear flowFor PPVF with non-Darcy flow effect, the fracture behaves like FPVF with the same height (hf) during this period, so the solution is corrected as follows:
(23)hfDpwDtD=2πkmhfpi−pwqrefμ=πΓ5/42CfD,apptD1/4,where CfD,app is a correlation for the apparent fracture conductivity, which is the function of initial conductivity (CfDi), the Reynolds number (FND), and the minipermeability plateau (kmr). After simulating numerous cases, the relationship between apparent and initial (true) conductivities that is suitable for the special case using the parameters seen in Figure 5 is achieved, which is
(24)CfD,appCfDi=kmr+1−kmr1+2.51FND.


### 3.2.2. Formation-linear flow

During this period as shown in Figure6(a), formation fluid enters the fracture in a direction perpendicular to the fracture surface, but the pressure response within fracture is not considered again. Note that the flow pattern in the fracture exhibits line-shaped during bilinear and linear flow regimes, which is confirmed in Figure 7(a). Based on the analytical solution for FPVF presented by Gil et al. [41], both penetration ratio and non-Darcy effect are further taken into account in this work. Similar to the derivation of bilinear-flow solution, we recast analytical solution for formation-linear-flow period as follows:
(25)hfDpwDtD=πtD+αCfD,app,where α is a constant given by the following: α=π/3 for CfD≥25, α=0.944 for 10≤CfD<25, and α=0.902 for 5≤CfD<10 [41].Figure 7
Pressure fields within two-dimensional fracture at (a)tD=10−6 and (b) tD=10−2.
(a)(b)


### 3.2.3. Early-radial (ER) flow

If penetration ratio is small enough,hfD should be smaller than 0.5 suggested by Al-Rbeawi and Tiab [11]; an ER flow regime is developed around and along the fracture length in YZ plane. The ER flow regime is identified by a horizontal line on pressure derivative plot with
(26)tD×pwD′ER=0.5.During this period, fracture behaves like a finite-conductivity horizontal well with axis along the formation boundary, and the reservoir fluids flow radially inYZ plane towards fracture as seen in Figure 6(b). The effective penetration length becomes fracture length rather than formation thickness or fracture height, which is given by
(27)pwD,ER=pwD2hD=2πkmLfpi−pwqrefμ.The corresponding horizontal line on derivative curve is corrected as follows:(28)tD×pwD,ER′ER=0.52hD.Note that the flow pattern deviates from the line shaped, becoming the radial shaped as shown in Figure7(b).


### 3.2.4. Compound-linear (CL) flow

CL flow often occurs after ER flow when both upper and bottom impermeable boundaries are felt. The fluid flow is predominantly perpendicular to the fracture face inXY plane, and the response does not reach the region beyond fracture tip. The PFVF behaves like a FFVF during this period, as shown in Figure 6(c). The characteristic of this flow behavior is very similar to the compound linear flow for multistage-fractured horizontal well presented by Chen and Raghavan [4] and fractured inclined well presented by Dinh and Tiab [10], which is identified by an approximated 1/2-slope straight line on log-log plot of pressure derivative.


### 3.2.5. Pseudoradial (PR) flow

During this period, the reservoir flow begins to converge radially to the partial fracture. The flow inXY plane dominates compared with YZ plane; alternatively, the flow in z-direction can be ignored. The PR flow is characterized with tD×p’wD=0.5 on the log-log plot.Depending on the reservoir, fracture, and non-Darcy flow parameters used in a single test, the duration of some flow regime might elongate or shorten, even disappear. The non-Darcy flow effect leads to extra pressure drop, but cannot alter the characterization of flow regimes identified by the Darcy behavior.


## 3.2.1. Bilinear flow

The fluids beyond the fracture face in vertical direction are not developed as show in Figure6(a). An explicit solution for a FPVF was provided by Cinco-Ley and Samaniego [40] during bilinear flow:
(22)pwDtD=πΓ5/42CfDitD1/4.Figure 6
Schematic of some special flow regimes occurring during the production for a PPVF.
(a)
Bilinear and linear flow(b)
Early-radial flow(c)
Compound-linear flowFor PPVF with non-Darcy flow effect, the fracture behaves like FPVF with the same height (hf) during this period, so the solution is corrected as follows:
(23)hfDpwDtD=2πkmhfpi−pwqrefμ=πΓ5/42CfD,apptD1/4,where CfD,app is a correlation for the apparent fracture conductivity, which is the function of initial conductivity (CfDi), the Reynolds number (FND), and the minipermeability plateau (kmr). After simulating numerous cases, the relationship between apparent and initial (true) conductivities that is suitable for the special case using the parameters seen in Figure 5 is achieved, which is
(24)CfD,appCfDi=kmr+1−kmr1+2.51FND.


## 3.2.2. Formation-linear flow

During this period as shown in Figure6(a), formation fluid enters the fracture in a direction perpendicular to the fracture surface, but the pressure response within fracture is not considered again. Note that the flow pattern in the fracture exhibits line-shaped during bilinear and linear flow regimes, which is confirmed in Figure 7(a). Based on the analytical solution for FPVF presented by Gil et al. [41], both penetration ratio and non-Darcy effect are further taken into account in this work. Similar to the derivation of bilinear-flow solution, we recast analytical solution for formation-linear-flow period as follows:
(25)hfDpwDtD=πtD+αCfD,app,where α is a constant given by the following: α=π/3 for CfD≥25, α=0.944 for 10≤CfD<25, and α=0.902 for 5≤CfD<10 [41].Figure 7
Pressure fields within two-dimensional fracture at (a)tD=10−6 and (b) tD=10−2.
(a)(b)


## 3.2.3. Early-radial (ER) flow

If penetration ratio is small enough,hfD should be smaller than 0.5 suggested by Al-Rbeawi and Tiab [11]; an ER flow regime is developed around and along the fracture length in YZ plane. The ER flow regime is identified by a horizontal line on pressure derivative plot with
(26)tD×pwD′ER=0.5.During this period, fracture behaves like a finite-conductivity horizontal well with axis along the formation boundary, and the reservoir fluids flow radially inYZ plane towards fracture as seen in Figure 6(b). The effective penetration length becomes fracture length rather than formation thickness or fracture height, which is given by
(27)pwD,ER=pwD2hD=2πkmLfpi−pwqrefμ.The corresponding horizontal line on derivative curve is corrected as follows:(28)tD×pwD,ER′ER=0.52hD.Note that the flow pattern deviates from the line shaped, becoming the radial shaped as shown in Figure7(b).


## 3.2.4. Compound-linear (CL) flow

CL flow often occurs after ER flow when both upper and bottom impermeable boundaries are felt. The fluid flow is predominantly perpendicular to the fracture face inXY plane, and the response does not reach the region beyond fracture tip. The PFVF behaves like a FFVF during this period, as shown in Figure 6(c). The characteristic of this flow behavior is very similar to the compound linear flow for multistage-fractured horizontal well presented by Chen and Raghavan [4] and fractured inclined well presented by Dinh and Tiab [10], which is identified by an approximated 1/2-slope straight line on log-log plot of pressure derivative.


## 3.2.5. Pseudoradial (PR) flow

During this period, the reservoir flow begins to converge radially to the partial fracture. The flow inXY plane dominates compared with YZ plane; alternatively, the flow in z-direction can be ignored. The PR flow is characterized with tD×p’wD=0.5 on the log-log plot.Depending on the reservoir, fracture, and non-Darcy flow parameters used in a single test, the duration of some flow regime might elongate or shorten, even disappear. The non-Darcy flow effect leads to extra pressure drop, but cannot alter the characterization of flow regimes identified by the Darcy behavior.


## 3.3. Sensitivity Analysis

According to Eq. (2), the effect of non-Darcy flow is determined by two dimensionless parameters (i.e., kmr and FND) with the Barree-Conway model. It is assumed that the formation thickness, fracture length, and initial conductivity are given, i.e., hD=1, hfD=0.1, and CfDi=100π. To highlight the improvement of this work, the effect of dimensionless parameters with the Barree-Conway model and the penetration ratio of fracture are mainly discussed here because effects of dimensionless thickness, fracture conductivity, fracture location, and asymmetrical wellbore have been extensively investigated and analyzed [17, 25].Figure8 shows the influence of non-Darcy flow effect on transient pressure behavior of PFVF. Figure 8(a) displays the pressure and derivative curves for FND=0.1, 1, 10, and 100 at kmr=0.1. The Reynolds number (FND) is proportional to production rate according to Eq. (3). Increasing FND represents the increase of production rate. Besides, the case that FND=0 or kmr=1 indicates Darcy’s flow equation. The limiting Darcy cases are introduced, which are denoted by blue dashed lines. Case 5 indicates the pressure responses for the initial conductivity in the Darcy flow condition, while case 6 indicates the minimum conductivity in the Darcy flow condition (CfDmin=CfDi×kmr). For case 1~4, non-Darcy effect causes a significant reduction in fracture conductivity in simulated cases. As a result, an extra pressure drop is caused by non-Darcy effect. Case 4 exhibits a whole sequence of flow regimes similar to the initial case (case 5) due to a weak non-Darcy effect. The increasing of FND shortens the duration of formation-linear flow regime until disappears (i.e., case 1 and case 2). Meanwhile, the pressure drop and derivative curves also increase and approach the minimum case (case 6) because of strong non-Darcy effect. As a result, the curves of pressure responses with non-Darcy effect (i.e., case 1~case 4) are distributed within the range between the initial case (case 5) and the minimum case (case 6). Note that case 1 is almost overlapped with the minimum case for times of interest.Figure 8
Effect of the Reynolds number on (a) pressure response and (b) stabilized flux distribution.
(a)(b)The flux distribution stabilizes and remain unchanged when pseudoradial flow is established, which is shown in Figure8(b). Note that the flux distribution indicates the arithmetic average value of flux distribution among fracture panels. The integral with regard to xfD represents the volume of fluids entering the fracture. In case 6, most of the flow occurs from the half of fracture away from the wellbore. With the FND increasing, the characterization of fluid flow is approaching the unified case (case 6): the flow from the tips decreases, and more fluid comes from the half of fracture closer to the wellbore.Figure9(a) plots the curves of pressure responses for kmr=0.1, 0.2, 0.5, and 0.9 at FND=10. Different values of kmr represent different minimum cases under Darcy flow, which are denoted by blue dashed lines. The initial case is denoted by green dashed line. Solid lines indicate the cases with non-Darcy effect. The non-Darcy cases are located within the range between the initial case and the minimum case. As seen in Figure 9(a), as the value of kmr decreases, the pressure drop and its derivative would deviate from the initial case and approach the corresponding minimum case due to intense non-Darcy effect. When the value of kmr approaches unity, such as case 4, the pressure response of the minimum Darcy case is similar to the initial Darcy case. The apparent conductivity for case 4 is always in the order of magnitude for infinite conductivity during the process of conductivity degradation. As a result, the effect of non-Darcy flow tends to disappear. When the value of kmr is relatively smaller, the non-Darcy effect becomes significant. Figure 9(b) shows the stabilized average flux distribution along fracture. As the value of kmr decreases with a strong non-Darcy effect, the apparent conductivity decreases, and more fluid comes from the half of fracture away from the wellbore. This is different from the phenomenon described in Figure 8(b), because the minimum conductivity is different, not a unified value.Figure 9
Effect of the relative minimum permeability on (a) pressure response and (b) stabilized flux distribution.
(a)(b)As analyze from Figures8 and 9, the dimensionless parameters of larger FND and small kmr make the non-Darcy effect more significant. With the consideration of non-Darcy effect, the apparent fracture conductivity exhibits different conductivities in the Darcy flow condition when it actually has only one true value. Figure 10 further investigates the effect of penetration ratio in different conditions of initial conductivity. The value of FND and kmr is still the same as that in Figure 10, and the transient pressure and its derivative show a similar trend. For smaller conductivity (CfDi=10π) as shown in Figure 10(a), the pressure drop is always lower than the minimum Darcy case but higher than the initial Darcy case all the time in the condition of hfD=0.1. However, the non-Darcy case would overlap with the initial Darcy case in the late-time period in the condition of hfD=1. As a comparison, a higher conductivity is considered as shown in Figure 10(b). The pressure response of non-Darcy case deviates from the initial Darcy case, but would overlap with the initial case on the onset of ER flow regime in the condition of hfD=0.1. The pressure response of non-Darcy case might overlap with the initial case on advance in the condition of hfD=1.Figure 10
Effect of the penetration ratio on pressure response in the condition of (a) low and (b) high conductivity.
(a)(b)In summary, the effect of non-Darcy flow has a more significant influence on the condition of smaller penetration and lower initial conductivity. Alternatively, the effect of non-Darcy flow has a more significant influence on the condition of small penetration.


## 3.4. Application to Field Case

The case presents a pressure drawdown test in an oil well in the Tahe oil field of Xinjiang, China. The vertical well is hydraulically fractured, and the fracture is recognized as a partial penetrating vertical fracture according to the results of microseismic monitoring. The basic parameters used in this case are summarized in Table1.Table 1
Reservoir and fluid properties for the well.
ParameterValueUnitFormation thickness,h32.5mVolume factor,B1.05Formation porosity,φ19.8Fluid viscosity,μ2.0mPa•sWellbore radius,rw0.1mTotal compressibility,ct4.7×10−4MPa-1Production rate,qw10.0m3/dFigure11 provides the pressure drawdown and its derivative data on log-log plot. According to the flow regime identification in Section 3.2, it is shown from the derivative curve that a wellbore-storage effect is firstly recognized with a unity slope. As a result, the bilinear flow regime is masked by the wellbore-storage period. Then, the formation-linear flow regime appears, followed by a significant ER flow regime. The slope of ER flow is a negative value not zero, and the value of tD×p’wD is approximately 0.25. According to Eq. (28), the penetration ratio is about 1/4, so the fracture height is about 8.125 m height. Meanwhile, the negative slope for ER flow regime satisfies the characteristic of low and moderate conductivity of PPVF as presented in Figure 10(a). Eventually, the compound-linear flow and pseudoradial flow are found successively.Figure 11
Pressure transient analysis for (a) flow regime identification and (b) type-curve matching of the pressure drawdown data.
(a)(b)To match the entire drawdown data, the effects of wellbore storage capacity and fracture damage should be incorporated into our model. Because we use Laplace-domain solution, the solution accounting for these effects can be flexibly written as follows:(29)p~owD=sp~wD+Sds+CDs2sp~wD+Sd.Through constantly calculating the unknown variables, Figure11(b) presents the measured and final simulated pressure drawdown data, and an excellent agreement exists between the measured and simulated pressure for drawdown test. This case uses the responses of several different flow regimes at the same flow rate, instead of the responses of single flow regimes at two different rates as presented by Guppy et al. [25]. Nevertheless, the non-Darcy flow responses can be still matched with the Darcy flow response for a lower conductivity because of the nonuniqueness problem. If non-Darcy effect is not considered in the matching denoted by dashed lines, the dimensionless conductivity (CfDi) is calculated as 125.56.If non-Darcy effect is considered denoted by solid lines, the dimensionless initial conductivity (CfDi) is calculated as 204.29 with kmr=0.187, and FND=71.56 all the time. The resulting matching parameters are given in Table 2. Considered that the production rate of this well is relatively high, the fracture might be subject to the effects of non-Darcy flow. The traditional analysis method would result in nearly 40% error in the estimation of conductivity. The type curves based on the Darcy model would yield an underestimated conductivity.Table 2
Type-curve matching results for the drawdown test.
Without non-Darcy effectWith non-Darcy effectParameterValueUnitParameterValueUnitWellbore storage,C0.0012m3/MPaWellbore storage,C0.0024m3/MPaSkin factor,Sd0Skin factor,Sd0Fracture conductivity,kfwf65.99D•mFracture conductivity,kfwf284.86D•mFracture length,Lf58.4mFracture length,Lf66.4mFormation permeability,km0.018DFormation permeability,km0.042D


## 4. Conclusions

Based on this work, a semianalytical solution was proposed for the pressure drawdown analysis of a well with partially penetrating vertical fracture. From the results of the investigation, several important conclusions need to be emphasized as follows:(1)
PPVF exhibited five typical flow regimes, including bilinear flow, formation-linear flow, early-radial flow, compound-linear flow, and pseudoradial flow. The corresponding analytical solutions can be presented by correlating with the penetration ratio and apparent conductivity(2)
The effect of non-Darcy flow makes the fracture conductivity behave lower than its true value. When the value ofkmr approaches unity or when the value of FND approaches zero, the effect of non-Darcy flow becomes weak, and the pressure transient responses approach the initial-Darcy case(3)
The non-Darcy flow has a more significant influence in the condition of smaller fracture penetration and lower initial conductivity. When the penetration ratio is smaller than 0.5, the pressure behavior exhibit a more remarkable variation with penetration ratio(4)
The type-curve matching of non-Darcy-flow responses with type curves without non-Darcy effect would yield lower estimate of true value of dimensionless fracture conductivity. As a result, fracture conductivity and fracture length would be underestimated



---

*Source: 8884750-2020-11-27.xml*