A mathematical model for estimating effective stimulated reservoir volume

This study presents a model application for the evaluation of Effective Stimulated Reservoir Volume (ESRV) in shale gas reservoirs. This current model is faster, cheaper, and readily available for estimating ESRV compared to previously published models. Key controlling parameters for efficient ESRV modeling, including geomechanical parameters and time, are considered for the model development. The model was validated for both single and multi-stage fractured reservoirs. For the single fractured reservoir, an ESRV of 3.07 × 106 ft3 was estimated against 3.99 × 106 ft3 of ESRV-FEM field data. Whereas, 7.00 × 109 ft3 ESRV was estimated from the multi-stage fractured reservoir against 7.90 × 109 ft3 of fractal-based model results. Stress dependence, time dependence, and permeability dependence of shale gas reservoirs are found to be essential parameters for the successful calculation of ESRV in reservoirs. An ESRV determined using this method can obtain the estimated ultimate recovery, propped volume, optimal fracture length, and spacing in fractured shale gas reservoirs.


Introduction
Hydraulic fracturing application in horizontal wells has significantly contributed to more production from shale gas reservoirs that were once difficult to exploit (Quainoo et al. 2020a). The source of this gas accumulation is the shale rock, and the geological trap is the shale rock itself (Guarnone et al. 2012). Shale gas production from these reservoirs requires a well-designed hydraulic fracturing treatment. Gas shales are believed to contain pre-existing natural fractures that interconnect with induced hydraulic fractures to form a complex fracture network known as stimulated reservoir volume (SRV) (Cipolla and Wallace 2014;Al-Rbeawi 2020). In terms of microseismic signatures, SRV is described as that volume of the shale reservoir which generates a high signal during fracturing treatment.
ESRV estimation is significant for the optimization of shale fields, design of fracture completion, and proper proppants placement. Moreover, an accurate study of the ESRV will lead to the design of a fracturing process that will stop unwanted damage to the shale formation (Jagaba et al. 2018;Abubakar et al. 2016). The existing methods for estimating SRV are; induced-fracture-monitoring, microseismic survey, and analytical or numerical modeling (Sheng et al. 2017;Hu et al. 2013). These methods, however, do not directly provide the volume of the reservoir that is producing rather, they give the volume of the reservoir that is affected by stimulation. Consequently, a new concept known as Effective Stimulated Reservoir Volume (ESRV) has been introduced to define the reservoir's volume responsible for production (Lacazette et al. 2015). Thus, the ESRV evaluation gives a clear indication of the productive volume. Numerous papers used varying nomenclature and methods in describing the ESRV. For example, some described it as Contributing Reservoir Volume (CRV) (Friedrich and Milliken 2013), Tributary Drainage Volume (TDV) (Doe et al. 2013), Drainage Volume (Rahman et al. 2020), and Effective Propped Volume (EPV) (Maxwell 2013).
Several methods have been used in the literature, such as the application of microseismic interpretation, fractal application, analytical estimation, and mathematical modeling, for estimating the ESRV. Microseismic interpretation is often used to obtain fracture distribution or geometry (Auger et al. 2011;Eisner and Stanek 2018). The fractal application means converting microseismic events into fractal geometry using Fractal Random-Fracture-Network Algorithm (FRFNA) (Sheng et al. 2017;Zhou et al. 2016). Also, integrating methods like the rate transient analysis (RTA) into the fractal application method have been implemented to study ESRV . Hydraulically fractured shale gas reservoirs can also be analytically deliberated to estimate the fractal ESRV (Zhang et al. 2018). However, mathematical models provide the premise for integrating geomechanical properties that drive the ESRV formulation (Zhang et al. 2018). Geomechanical models were some of the first applications used to evaluate the ESRV. Lately, mathematical models are continuously developed to compute the ESRV simply and cost-effectively ). In such models, the flow rate, fracture conductivity, and formation stress were the parameters used to obtain ESRV (Moos, et al. 2013). Downhole and surface seismic data have been used for careful data acquisition in this area. Moreover, the discrete fracture network (DFN) modeling was applied to simulate the ESRV for either the whole fractured area or individual fracture stages (Doe et al. 2013;Danso, et al. 2021). It monitors the parameters for gas flow from each distinct fracture. These applications have found that the key controlling parameters involved are the formation stress, hydraulic fracture propagation, reservoir pressure, fluid flow, and type of rock failure (Quainoo et al. 2020b). So, modeling these processes using previous techniques is complex as it involves incorporating the various modules of each parameter. In some of the methods, the geomechanical parameters and time dependence of the production are neglected. Besides, microseismic data acquisition and interpretation are often costly and time-consuming (Eisner and Stanek 2018;Hussain et al. 2017).
Therefore, this paper intends to present a simple model for estimating the ESRV in both single-stage and multistage fractured shale gas reservoirs. The vital geomechanical variable such as stress (stress dependence), bulk modulus, and poisson ratio are considered in this model. A simulated hydraulic fracture geometry that can be used to estimate the ESRV is also introduced. Table 1 shows a review of the literature detailing the contributions and limitations of various ESRV methods.

Proposed mathematical model
The mathematical model consists of three major components in the model which are constituted from the fracture length, fracture height, and geomechanical variables (stress, poisson ratio, bulk modulus, etc.). They are described as: • Effective fracture length • Effective fracture height • Distance of Investigation Moreover, the distance of investigation is considered in the mathematical formulation to define the extent at which the reservoir is contributing to gas production. Alternatively, the fracture geometry is represented by both the effective fracture length and effective fracture height. These represent the parts of the fracture geometry that contributes to production. The hydraulic fracture propagation can be estimated with respect to time considering the effective fractures. Each component is carefully selected and described to finally arrived at the new approach and simple mathematical model for ESRV estimation.
The validation was performed by considering two case studies as explained in chapter 2. These case studies are that of: • A Fractal Random-Fracture-Network algorithm (FRFNA) and • An ESRV fractal evaluation model (ESRV-FEM).

Case 1: the fractal fracture network algorithm-based model
The new ESRV model would be validated against an FRFNA-based model that used microseismic data to determine the ESRV (Sheng et al. 2017). In their model, a Fractal Random-Fracture-Network algorithm (FRFNA) was employed to convert Microseismic data into fractal geometry. Also, an orthogonal fracture coupled dual-porositymedia representation elementary volume (REV) flow model was introduced to foresee the volumetric flux in shale reservoirs. The ESRV was determined by an ESRV criterion, which assumed the SRV region's random distribution in both porous kerogen and inorganic matrix. The FRFNA-based model also attempted to simplify the complex fracture by converting microseismic results into regular pattern fractal geometry. The fractal fracture network nodes are matched in order to calibrate the fracture geometry using the integer programming. It has 12 steps to estimate the ESRV.
(1) Upload microseismic data and extract the coordinate in MATLAB and (2) judge the bifurcation direction base on the distribution features of fractures density (forward or reverse bifurcation). In this step, forward-bifurcation tree-generating rules or reverse-bifurcation tree-generating rules will be determined, (3) propose iteration times to control the extension of the fractal fracture. Initial iteration times the value of 2 is assigned, (4) An approximate range of generating length according to the distribution of microseismic data is proposed. The range is divided into 100 parts. Initial loop variable for generating length (i) to 1 is set, (5) An approximate range of deviation angle according to the distribution of microseismic data is proposed. The range is divided into 100 parts, and the initial loop variable for deviation angle (j) to 1 is set, (6) set loop the variable for generating rules (k) to 1, (7) introduce the random function to determine specific generating rules, (8) the corresponding fractal geometry by combining times id generated, generating length, deviation angle, and specific generating rules, (9) the distance between every microseismic datum and the nearest fracture node, and the distance between every microseismic datum and the nearest fracture heel is summed. If the matching rate is larger than the given criterion (90% in the paper). When the matching rate is smaller than the given criterion, if loop variable for generating rules (k) is smaller than 100, then add one to k and repeat step. Otherwise, move on to the next step, (10) one to the loop variable for the deviation angle (j) is added. If j is smaller than 100, repeat step 5. If not, move on to the next step, (11) add one to the loop variable for the generating length (j). If i is smaller than 100, repeat step 4, if not, move on to the next step and (12) 1 is added for initial iteration times, repeat step 3.

Case 2: ESRV-fractal evaluation model
The multi-stage fracture reservoir, on the other hand, was validated using ESRV-FEM model reported by Zhang et al. (2018). In the ESRV-FEM model, multiple transport mechanisms, SRV, and unstimulated SRV were incorporated in the model. The heterogeneity extent identified in fractal indexes were used to compute the ESRV at various coefficients of inter-porosity flow and storage ratio. As permeability and porosity distributions are complex, the ESRV-FEM recognizes a simplified form of the fracture. It distinguishes three flow regions such as hydraulic fracture, stimulated reservoir volume, and unstimulated reservoir volume. It also assumed a constant initial pressure and gas production rate. Moreover, it adds up diffusion equations with initial and boundary conditions. So, the SRV line source is derived and integrated along the hydraulic fracture to obtain the flowing bottom hole pressure. This corroborate with modeling the extend at which gas is produced from the reservoir in the methodology section.
These studies were selected as they considered the fracture geometry pattern to arrive at the ESRV. They are also the latest in the field. In addition, the difficulty in data acquisition is made easier using this model as well as validating it. However, the main ESRV controlling parameters are involved in the formulation for accurate computation. The data used for the FRFNA-based model is presented in Table 2. It is from the Sichuan Basin of the fuling field, China. Likewise, the data for ESRV-FEM validation are from the unconventional gas reservoirs in China. It is shown in Fig. 3.
Different papers have prominently presented some procedures to calculate ESRV from their available field data indirectly. However, they did not explicitly provide a direct method of determining the ESRV. As a result, this methodology intends to offer an easy-to-use mathematical model that can estimate the ESRV. For the mathematical model development, the ESRV is intended to be evaluated by considering the effective fracture length, effective fracture height, and distance of investigation (DOI) as in Fig. 1. The multiplication of fracture length, fracture height, and width makes up the fracture geometry that promotes the passage of the shale gas. The effective fracture length (l eff ) and height (h eff ) are the productive paths after permeability and stress impact on the fractured shale reservoir. They provide the conduits that are enough for gas to escape from the shale gas reservoir. The current model is expected to be faster, cheaper, and reliable as it considers the treatment process's stress and time dependency. In this model, therefore, geomechanical parameters such as stress, poisson's ratio, and modulus of elasticity are considered.
Furthermore, this work intends to visualize the fracture geometry using the pseudo-3D model in MATLAB software to determine the effective fracture width and height. The P3D model is compared with the well-known KGD and PKN models. This fracture propagation model is used to illustrate the fracture geometry (the width, w f and height, H f ). The fracturing process is anticipated to propagate perpendicular to the minimum horizontal stress, Shmin. A cartesian coordinate system is established to show the fracture as in Fig. 1a below: This method of finding ESRV is divided into two subsections: ESRV in single-stage fractured reservoirs and multistage fractured reservoirs. For a single-stage fractured reservoir, ESRV is calculated using Eq. (1).
ESRV for multi-stage fractured reservoir can be calculated using Eq. (2).
The fracture geometry is a complex function of initial reservoir conditions, heterogeneous reservoir, and geomechanical properties, permeability, porosity, natural fractures, and treatment factors such as pressure, volume, and injection rate. The above-simplified assumptions were made to simplify this complex scenario by depending on fracture geometry's influence in defining fractured shale gas reservoir. So, the well perforations can be delineated by the extent of the hydraulic fracture propagation. The well performance can be measured directly through this fracture geometry-based approach. Therefore, this model will use the extent of the fracture geometry to define the ESRV. The main controllable factors that characterize the geometry are: injected fluid volume, fluid loss, injected fluid viscosity, and injection rate. Fracture width can be increased by increasing the viscosity of fracturing fluid and injection rate for an optimized penetration distance. This work provides a direct perspective on how to improve the ESRV for better-fractured shale gas reservoir performance.

Effective fracture length
Fracture length is the most crucial because it affects the effectiveness of the fracturing process (Yu et al. 2014). Fracture length is the distance from one tip of a fracture to another. Often, the term 'half-length' is used to represent the distance between the well and tip of the fracture and half of the length of the fracture. The fracture length has a significant impact on the hydrocarbon flow rate and cumulative production. Both the flow rate and the production increase with an increase in the fracture length. However, longer fracture increases the operational cost due to the need for more proppants, additives, etc. An optimum fracture length is therefore required for economic hydrocarbon production (Zhang et al. 2009).
Effective fracture length entails the whole length of the fracture that contributes to the production of gas. It is obtained by applying the explicit expression of the fracture length in Eq. 3. It is suitable for this model as it assumes equal flow rate distribution in the fractures and constant net pressure. The expression is given as follows (Wu 2018): where Q constant injection rate (ft 3 /h). E′ Elasticity plain strain modulus (psi).
Moreover, constant net pressure and equal flow rate distribution assumptions were made along the fracture which is applied to express the fracture length and width formula. The hydraulic fracture length and fracture width at the wellbore are simulated to portray their evolution in the reservoir as in Figs. 6 and 7.

Effective fracture height
Fracture height is one of the parameters that govern the geometry of any hydraulic fracturing. It is evaluated using the net pressure of fracture (P net ) and fracture toughness K IC . In hydraulic fracturing, the loading conditions and fracture geometry are assumed to be simple, as there is inadequate knowledge about the fracture shape peculiarities. Because the fracture propagates according to linear elastic fracture mechanics, the stress field near the crack tip is evaluated using the elastic theory. Thus, if the stresses near the crack tip exceed the toughness of the fracture, the crack grows (Valko and Economides 1995). In this improved model, it is crucial to consider the geomechanical properties in determining the effective fracture height. With the net pressure and fracture toughness, the effective fracture height can be evaluated using the fracture height expression as follows (Ren et al. 2018): where K IC fracture toughness (Pa m 1/2 ). P f Fracture fluid pressure (Pa). min minimum horizontal stress (Pa).

Distance of investigation (DOI)
When the fractured reservoir continues production at a constant rate, the wellbore pressure continues to diminish as time increases. Simultaneously, the size of the drained reservoir area increases, and the pressure transient adjusts farther out into the reservoir. The effective distance by pressure transient traveled is referred to as the DOI (Behmanesh et al. 2015). In other words, it represents the depths at which the properties of the formation are being examined during a test at any time. Therefore, DOI is the distance at which the pressure reduction is negligible.
A very comprehensive analysis of the DOI for radial geometry was performed by Kuchuk (2009). After that, many expressions have been reported for determining DOI (Zheng 2016). However, none of them has been used for linear flow except that of the equation reported by Behmanesh et al. (2015) given in Eq. (3). Thus, the equation by Behmanesh et al. (2015) was used in the current model.
w h e r e C 1 . 4 1 ( f o r l i n e a r f l o w ) , k = Formation permeability (md) , = Porosity (% ) , c t = Total compressibility psi −1 , g = Gas viscosity (cp) , t = time (h).

Model validation and discussion
The single-stage fracture model was validated against a FRFNA-based model results. It used microseismic data to determine the ESRV (Sheng et al. 2017). In their model, a Fractal Random-Fracture-Network algorithm (FRFNA) was applied to transform Microseismic data into fractal geometry. Moreover, an orthogonal fracture coupled dual-porositymedia representation elementary volume (REV) flow model was presented to provide the volumetric flux in shale gas reservoirs. The ESRV was defined by an ESRV criterion, which presumed the SRV region's random distribution in both porous kerogen and inorganic matrix. On the other hand, the multi-stage fracture reservoir was validated using ESRV-FEM model results reported by Zhang et al. (2018). In the ESRV-FEM model, SRV, multiple transport mechanisms, and unstimulated SRV were incorporated in the model. The heterogeneity extent identified in fractal indexes was used to compute the ESRV at various inter-porosity coefficients with a storage ratio. The field data used for FRFNA-based and ESRV-FEM models are listed in Tables 2 and 3, respectively.
The results from the proposed model are compared with the previous model's results as shown in Table 4. The estimated fracture length is plotted in Fig. 2, and it follows the same trend as the analytical solutions (KGD and PKN models) reported by Zadeh (2014) as it propagates over time into the reservoir. From Table 4, the improved model resulted in an ESRV value with a decrease of at least 12%. This supports the assertion that microseismic data tend to overestimate the productive volume in shale gas reservoirs. In addition, the current model has included the basic and some of the most important geomechanical parameters in both fracture height and fracture length modules in order to obtain a more reliable ESRV estimation of both reservoirs. Besides, the stress and permeability sensitivity were considered in this model. Thus, this model can be regarded as acceptable to be used for ESRV estimation of shale gas reservoirs since it gave values that are close to those reported in the published pieces of literature.
In Fig. 3, the time dependence of the ESRV is illustrated. The fracturing process is fast and time-dependent. Other models assume constant time for all the fracturing treatment regardless of the formation's geomechanics (Valko and Economides 1995). This proposed model considers the geomechanics and therefore introduces a time variable. Our result gives the opportunity for any fracturing process to establish the ESRV as fast as possible concerning time which will assist in optimizing the treatment. The time shows the variation of the ESRV during fracturing, which other models do not provide. The model can be applied easily to show the extent of the productive volume and can be used to optimize the stimulation process before or during fracturing. As a single-stage fracture, it provides the expected values that made up the effective fracture geometry. This geometry controls most of the gas production in shale gas reservoirs. In addition, the ESRV recorded by applying the improved model is lower than the FRFNA-based model as illustrated in Fig. 3. The FRFNA-based model computed the ESRV using fractal geometry that is derived from microseismic data. This data generally records parts of the reservoir that are not connected to the fractures that are responsible for production. So, the FRFNA-based model overestimates the ESRV. In other words, the model depends on microseismic data to provide the fracture distribution, which tends to overestimate the extent of the ESRV. Whereas the improved model incorporated the significant geomechanical variables (such as stress, Poisson's ratio, bulk modulus, etc.) that should describe the productive volume of a fractured shale gas reservoir. Therefore, these variables provide better result that delineate only the productive volume of the  Fig. 3 Comparison between the proposed model for the single-stage fracture model and the FRFNA-based model fractured reservoir. The complex fractures generated can be easily optimized while considering these variables as they help to provide the realistic ESRV and control the stability and/or deformation of the shale formation. Also, as stated earlier, these are part of the essential parameters that control gas production in fractured shale gas reservoirs. Furthermore, from Table 2, the ESRV-FEM, which was made up of multiple transport variable, SRV and unstimulated SRV, provides lower ESRV than the improved model. The data used is for a multi-stage fractured reservoir. Hence, the cumulative effective fracture lengths and heights are larger than that of single fractured reservoirs. The ESRV obtained from the improved model was 7.0 Bcf which is quite close to 7.9 Bcf obtained with ESRV-FEM. The reason for this increase could be as a result of the incorporation of the unstimulated SRV in the ESRV-FEM model. From both Figs. 4 and 5, most of the gas will be produced at an early stage from the ESRV. In order to sustain high gas production from the ESRV, optimal fracture spacing is needed during the fracture treatment. The effect of stress shadow in multi-stage fracturing may reduce gas production gradually or sharply with time. It depends on how well the geomechanical properties are considered in the estimating model. So, this model considered the stress dependence of the reservoir in its formulation and assumptions. Besides, a hydraulic fracture simulation of the single fracturing stage was presented, and it corresponds with the model's fracture geometry result. It is believed that the result obtained shows how efficient it can be in evaluating ESRV in shale gas reservoirs (Figs. 6, 7).

Conclusion
An improved model is proposed to estimate the effective stimulated reservoir volume (ESRV) in shale gas reservoirs through the fracture geometry variables. It considers the stress and time dependences of the fractured shale gas reservoirs. The improved model was validated in single-stage and multi-stage fractured reservoir. The proposed model shows how ESRV varies with time due to stress and/or permeability changes. So, the model can be used to predict the ESRV over a period. A decrease in volume from improved models confirms that microseismic events overestimate the productive volume. Thus, incorporating the geomechanical parameters gives a more accurate model of the volume of the shale gas reservoir that is productive. It is relevant to include shale gas reservoirs' stress and time dependencies in any model built for the purpose of ESRV computation in such reservoirs. Subsequently, more analysis can be made with this model on any fractured shale gas reservoir to define the ultimate recovery estimate, optimal fracture length, fracture spacing, and propped volume of the stimulated reservoir. Moreover, a laboratory and/or numerical simulation may be applied to validate the improved model.