Abstract
Hydraulic fracturing enhances hydrocarbon production from low permeability reservoirs. Laboratory tests and direct field measurements do a decent job of predicting the response of the system but are expensive and not easily accessible, thus increasing the need for robust deterministic and numerical solutions. The reliability of these mathematical models hinges on the uncertainties in the input parameters because uncertainty propagates to the output solution resulting in incorrect interpretations. Here, I build a framework for uncertainty quantification for a 1D fracture geometry using Woodford shale data. The proposed framework uses Monte-Carlo-based statistical methods and is comprised of two parts: sensitivity analysis and the probability density functions. Results reveal the transient nature of the sensitivity analysis, showing that Young’s modulus controls the initial pore pressure, which after 1 hour depends on the hydraulic conductivity. Results also show that the leak-off is most sensitive to permeability and thermal expansion coefficient of the rock and that temperature evolution primarily depends on thermal conductivity and the overall heat capacity. Furthermore, the model shows that Young’s modulus controls the initial fracture width, which after 1 hour of injection depends on the thermal expansion coefficient. Finally, the probability density curve of the transient fracture width displays the range of possible fracture aperture and adequate proppant size. The good agreement between the statistical model and field observations shows that the probability density curve can provide a reliable insight into the optimal proppant size.
1. Introduction
Hydraulic fracturing (HF) improves the flow rate of recoverable reserves from low permeability geological formations [1, 2]. The process induces rock failure by simply injecting fluid into the reservoir to counteract the in situ stresses and its tensile strength. Such newly created fractures are maintained using proppants often mixed with the fracturing fluid. HF is an expensive practice that requires careful design and good numerical modelling, and a principal problem with these mathematical models is that they are based on ideal systems with physical simplifications, thus neglecting either fracture initiation from wellbore, or near-tip effects, or fracture aperture evolution of preexisting fractures [3]. Although these simplifications facilitate the problem-solving process, they introduce sources of error that reduce the degree of confidence in the output solution [4]. Furthermore, uncertainty is embedded in the initial conditions, boundary conditions, geometry, and input parameters [4]. Some uncertainties create a large variance of the output solution and contribute most to the prediction of the degree of confidence while others have smaller effects [5]. Here, I address the uncertainties in shale formation properties generally calculated in a laboratory/field setting. Shale is a sensitive material, easily disturbed by normal drilling techniques, and specimens are difficult to sample, core, characterize, and test [6]. Furthermore, well logs measure formation properties such as Young’s modulus, Poisson’s ratio, and permeability, but well log interpretation is a subjective practice, where the tools have a margin of error and the measured properties are spatially limited. In this study, I investigate Woodford shale, a major source rock located in the Midwestern U.S. Its low permeability makes it difficult to fully profit from its potential, making it a perfect case study for HF. Field observations for Woodford shale [7] show that its natural fractures have a vertical dip and strike West-East with an average distance between fractures of 1.2-1.5 (m) (Figure 1). Hydraulic fracturing reactivates these preexisting fractures, thus creating multibranch fractures (Figure 2).
Quantifying the effects of spatial variability in formation properties [8, 9] on the reliability of hydraulic fracture simulations has been studied [6, 10–13] but is restricted by simplified deterministic solutions or computational timescales of numerical solutions. The uncertainty quantification for the simple linear elastic model given by [12] calculates the range of possible fracture apertures using Monte-Carlo simulations where Young’s modulus and the confining stress are random parameters, and they consequently reformulate the governing equations as stochastic partial differential equations (SPDEs). They adopt a stochastic collocation method to solve these SPDEs and replace the classic Monte-Carlo methods with two different methodologies: (1) the input parameters respect a log-normal distribution function, where only the mean and the dispersion values are used to calibrate the model, and (2) the input parameters fall within 96% of the log-normal distribution (a method called maximum entropy). They show that the probability distribution function of Young’s modulus and the confining stress greatly impact that of the fracture aperture. Furthermore, [12] were inspired by [8] who developed the stochastic framework to efficiently perform uncertainty quantification for fluid transport in porous media in the presence of both stochastic permeability and multiple scales and present the solution as a polynomial approximation. In a recent publication, [13] address variability in discrete fracture network geometry using the very robust multilevel Monte-Carlo methods, rather than the stochastic methodology claiming that the classic stochastic collocation fails to provide reliable estimates of first-order moments of the output solution due to the lack of smoothness. In another recent publication, [14] proposes the (MLMC) methods in conjunction with a well-assessed underlying solver to perform DFN flow simulations using Darcy flow and prove that the method is robust enough to tackle complex geometrical configurations, even with very coarse meshes. [15] address the heterogeneities in fractured reservoirs and assume that the fracture intensity, orientation, and conductivity of different fracture sets are the uncertain parameters. They confer each parameter with a probability distribution and create an integrated history matching workflow for fractured reservoirs assuming numerical methods for multiphase flow simulation, while updating the fracture properties via dynamic flow responses, such as continuous rate and pressure measurements. They generate initial realizations using Monte-Carlo simulations and based on observed fractures at the well locations. The automated history matching approach results in multiple equally probable discrete fracture network models, where upscaled flow simulation models honour both the geological information and the dynamic production history. Furthermore, [16] use uncertainty quantification for coupled subsurface flow and deformation processes. They circumvent the slow convergence of the traditional Monte-Carlo methods using an intrusive polynomial chaos expansion method for Biot’s poroelasticity equations based on Galerkin projection with uniform and log-normally distributed material parameters. Results show good agreement with the Monte-Carlo and ANOVA-based probabilistic collocation methods.
Although significant studies have been presented on uncertainty quantification in hydraulic fracturing, they are restricted to finite element schemes with defined mesh size and number of random values. Furthermore, they focus on geometric properties or are limited to elastic formation parameters, rather than rank all formation properties including the thermal properties, which is covered in this paper. The purpose here is to address formation inhomogeneity and spatial variability in low permeability shale reservoirs and their impact on fracture response using a unique but simple framework for uncertainty quantification comprising of sensitivity ranking and probability density curves. Then, I test the validity of the framework on the 1D porothermoelastic analytical solution for pressure, temperature, leak-off, and width given by [3], with application to Woodford shale data. A series of simulations demonstrate uncertainty propagation in the deterministic solution and describe the probability density function (PDF) for the fracture width to explain why some proppants have higher success rates than others.
2. Materials and Methods
2.1. Conceptual Model
Analytical solutions for pressure, temperature, strain, and displacement are derived in detail by [3] following the 1D consolidation theory.
Figure 3 zooms in to a multibranch fracture and shows how the fracture apertures are perpendicular to the minimum horizontal stress and intersect a porous medium characterized with an initial pore pressure and an original temperature and subject to fracturing fluid with pressure and temperature at its faces (Figure 3). Figure 3 also shows that the distance between two consecutive fracture apertures is , but because of the mirror symmetry, the one-dimensional solution is constrained to the length . Furthermore, the pore pressure gradient induces fluid leak-off from the fracture into the formation at each fracture face, and impermeable, adiabatic, and frictionless boundary conditions prevail at the axis of symmetry.
2.2. Governing Equations
The one-dimensional consolidation for transversely isotropic material is described in detail by [3], but a brief summary of the equations is presented here for the sake of completeness. Assuming the () plane is the isotropic plane and the axis is the axis of rotational symmetry (Figure 4), the coupled equations for pressure, temperature, and displacement represent the response of the hydraulic fracture system.
First, the analytical pressure solution is derived as then, the temperature is given by finally, the displacement is given by where is the half distance (m), is (m-1), is the stress boundary condition that is the difference between the fracturing fluid pressure (MPa) and the minimum horizontal stress (MPa), is the stiffness coefficient (MPa) calculated from Young’s modulus and Poisson’s ratio, is a thermal parameter calculated from the stiffness coefficients and the solid thermal expansion coefficient () and is Biot’s coefficient [1], and , , and are lumped parameters defined as and , and are defined as where is the fluid viscosity (Pa s), is the permeability (m2), is a poroelastic parameter in (MPa), is the thermal conductivity , and is the overall thermal capacity (). Note that subscripts 1 and 3 represent the direction and the direction, respectively.
The initial and boundary conditions associated with Equations (1)–(3) are given as follows:
Despite the complexity of the porothermoelastic solution, one can easily identify the contributing formation properties. The theory of linear elasticity assumes the rock as a dry solid material characterized by Young’s modulus (GPa), and Poisson’s ratio [1]. However, pore fluid alters the mechanical behaviour of the porous rock, thus introducing three additional poroelastic parameters: Biot’s coefficient [1], Biot’s modulus (GPa), and the mobility ratio (ratio between permeability of the formation and viscosity of the pore fluid: ) (). In addition, the one-dimensional anisotropic solution couples the poroelastic solution with the effects of the thermal gradient (between the fracturing fluid and the rock formation), thus introducing thermal parameters to the solution: the solid skeleton linear thermal expansion coefficient in the direction, the volumetric thermal expansion coefficient of the pore fluid , the thermal conductivity of the porous medium (in the direction) and the overall thermal capacity . Finally, the anisotropic nature of shale has a significant impact on its mechanical responses [3] and is included in this study. Table 1 enumerates the formation properties contributing to the uncertainty study for the 1D anisotropic porothermoelastic solution.
In the following, I detail the proposed uncertainty quantification (UQ) framework comprised of two parts: (1) sensitivity ranking to rank the input parameters by order of degree of influence on the output solution and (2) probability density function of the output solution that depicts all possible HF responses and their corresponding likelihood.
2.3. Part I: Sensitivity Analysis
To analyse the propagation of uncertainties, one must study how the output solution varies when the input variable takes on different values [5]. Using the popular Monte-Carlo methods, one generates a set of realizations of the input data , each characterized by a known probability distribution law, and calculates the correspondent solution of the model , then measures its statistics [4]. The Monte-Carlo methods have the advantage of being robust; that is, there are no constraints on the variances of the input data and the output. These methods are independent from the dimensionality of the random values (RVs) but require the deterministic solution and are characterized with a slow convergence rate. Although there are methods to increase the convergence rate such as variance reduction methods, a convergence rate of governs all Monte-Carlo methods, thus requiring a large population size [4, 5]. As a result, these methods are impeded by high computational cost and time.
Sensitivity analysis addresses the first part of the framework, where I rank input parameters by their order of influence on the variance of the solution. One of the simplest and most practical methods to perform sensitivity analysis is the one at a time method (OAT). In this method, one studies each parameter individually to minimize the chances of computer crashes. Typically, one random input parameter is assumed, while maintaining all others to their originally assigned values. The sensitivity analysis is driven by this random parameter that produces a host of possible output solutions. The steps for the OAT approach are listed in the following text for a deterministic solution, with input parameters. (1)Consider a random input parameter that fits a defined distribution function, e.g., normal distribution, with an assigned mean and standard deviation, and generate of . Figure 5 shows the population distribution for Young’s modulus (GPa) converging for (2)Calculate the deterministic solution for each value of (3)Calculate the variance of the values of output solution and note it as (4)Repeat procedure for each of the other input parameters (5)Calculate the sum of the output solution for values of variance. (6)Calculate the coefficient of variation for each input parameter. This coefficient represents the degree of influence of each input parameter on the output solution
2.4. Part II: Probability Density Curves
Proposed methods for uncertainty quantification include the probability density function (PDF), where its value is always positive and its integral over the entire range of possible values equals to one. To generate the PDF, one simultaneously introduces the high-ranking random input variables (from the first part of the framework) in the deterministic solution and fits the range of output results to a probabilistic distribution type (normal, triangular, Gumbel, lognormal, etc.). The procedure for a deterministic solution with input parameters is as follows: (1)Select the input parameters with a high sensitivity ranking (from sensitivity analysis) as random variables and keep the formation parameters with a negligible effect to their original value(2)Generate a population of values for each of the selected random variables and calculate the output solution(3)Plot the histogram of all the n output values (4)Fit the histogram to a distribution function and calculate the PDF function
3. Results
3.1. Woodford Shale Data
I test the proposed framework on the 1D solution with typical Woodford Shale data presented in Table 2 [17–20]. The transversely isotropic nature of this low permeability shale signifies that its properties depend on the direction, where has a subscript of 1 and a subscript of 3 (Figure 3), except for the linear thermal expansion coefficient that is independent of the direction and is replaced with one value . The model considers thermal effects because of the significant temperature difference between the fracturing fluid temperature and the formation temperature.
However, literature review [17, 19, 20] confirms a large variability in some of those properties, summarized in Table 3.
3.2. Sensitivity Analysis: Application to Woodford Shale
The saturated porous medium in Figure 4 is subjected to an instantaneous compression due to injection, and as a result experiences a pore pressure jump, also known as Skempton’s effect. The pressure at the shale-fracture interface is the constant injection pressure , and consequently, the sensitivity score for all parameters is zero at (m) (Figure 6). Figure 6 shows that the dominant control on the initial pore pressure (Equation (1)) is the elastic Young’s modulus with a sensitivity score of 90%. As pressure distribution stabilizes inside the half-space, the sensitivity of the pore pressure to decreases and other formation properties come into play, as shown in Figure 7. Figures 7(a)–7(c) show that as the pore pressure front diffuses inside the formation, the sensitivity of pressure to significantly decreases and the pore pressure becomes more sensitive to the mobility ratio , the thermal conductivity , and the skeleton linear thermal expansion coefficient , respectively. The reason why thermal properties influence the pore pressure response is that the temperature gradient forces the fluid and rock to contract as the cooling front penetrates further inside the half-space.
(a)
(b)
(c)
is the second most dominant parameter, where its influence reaches 0.09 (m) at (min) (Figure 7(a)), 0.12 (m) at (min) (Figure 7(b)), and 0.18 (m) at (min) (Figure 7(c)), and is the third most dominant parameter, where its influence reaches 0.06 (m) at (min) (Figure 7(a)), 0.09 (m) at (min) (Figure 7(b)), and 0.12 (m) at (min) (Figure 7(c)). The thermal conductivity sensitivity curve presents a peak that coincides with the dip in the mobility ratio sensitivity curve for (m) at (min) (Figure 7(a)), (m) at (min) (Figure 7(b)), and [m] at (min) (Figure 7(c)), highlighting a lower pressure zone, as seen in [3]. After 30 (min) of injection, the dip of the mobility ratio stabilizes at 48% and the peak in the thermal conductivity stabilizes at 38%. Figures 7(a)–7(c) also show that beyond the pore pressure front lies a zone where pore pressure is only sensitive to and , such as (m) at (min) (Figure 7(a)), (m) at (min) (Figure 7(b)), and (m at (min) (Figure 7(c)).
I proceed to perform sensitivity ranking for the temperature (Equation (2)), where contrary to the pore pressure, temperature is independent of any formation parameter at time , because the temperature front takes finite time to diffuse inside the formation.
Figure 8(a) shows that at (min), the cooling front reaches 0.11 (m) and depends on the thermal conductivity and the overall thermal capacity . Figures 8(b) and 8(c) show that the sensitivity of the temperature to increases while the influence of continues to decrease but remains dominant, with a sensitivity score of 70% for and for of 30% for at (min) and for (m). This is particularly interesting because while pressure depends on thermal properties, temperature only depends on two temperature properties, where the temperature at the shale-fracture interface ( (m)) only depends on the thermal conductivity Note that the Monte-Carlo simulation for temperature becomes smoother as the system stabilizes and reaches a steady condition.
(a)
(b)
(c)
Next, the 1D solution provides an expression for the fluid flow based on Darcy’s law, where the leak-off velocity is the flow at the shale-fracture interface ( (m)).
Figure 9 shows that the leak-off velocity is most sensitive to the mobility ratio and the thermal conductivity , with a sensitivity score of 90% and 10%, respectively.
Finally, the fracture width is the displacement at the shale-fracture interface, (m) (Equation (3)).
Similar to the pressure (Figure 6(a)), fracture width is initially () most sensitive to the elastic Young’s modulus (Figure 10), but because thermal effects are factored in the rock deformation process, the influence of decreases with time, while the linear thermal expansion coefficient the mobility ratio , and the thermal conductivity start to play a larger role in the displacement solution. Figure 10 shows that after 60 (min), the fracture width becomes most sensitive to (75%), followed by (12%), (9%), and (3%).
Table 4 summarizes the formation properties with a high sensitivity ranking. I consider the parameters with Y attribute for this second part and generate the probability density functions to quantify uncertainty in the output solution.
3.3. Probability Density Functions: Application to Woodford Shale
Figures 11(a)–11(c) show the probability density functions (PDFs) of the pore pressure response (considering the simultaneous variability of the four parameters indicated with Y in Table 4), at times (min), (min), and (min), respectively. The pore pressure variability fits a normal distribution function, where the pressure is fixed at (MPa) at the shale-fracture interface, (m). That is, PDF at this boundary is 100% and is not included in the plots.
(a)
(b)
(c)
Figure 11(a) shows that the transient pore pressure front diffuses from the boundary and reaches a depth of (m) at (min), to stabilize at a pressure range of 30-35 (MPa), for (m), which corresponds to the pore pressure jump felt in the porous medium at . Figure 11(b) shows that the pressure diffusion reaches 0.15 (m), and Figure 11(c) shows diffusion reaches 0.2 (m). Also, it is important to talk about the standard deviation because it indicates the amount of dispersion or variation observed in the population [21]. The pressure PDF follows a normal distribution function and has the highest standard deviation at (min), which decreases with time, where the PDF curve becomes sharper at 30 (min) and 60 (min).
Figures 12(a)–12(c) show the probability density functions (PDF) of the thermal diffusion, assuming a random thermal conductivity and a random overall heat capacity , calculated at times (min), (min), and (min), respectively. Like pressure, PDF of the temperature at the shale-fracture boundary is 100% (constant injection temperature of 30 (°C)) and is not included in the plots. The thermal diffusion front fits a normal distribution function and shows that the temperature of the porous medium decreases as a cooling front progresses deeper inside the medium to reach a depth of about 0.13 (m) at (min) (Figure 12(c)). The temperature also follows a normal distribution PDF, where its standard deviation, contrary to the pore pressure, does not seem to decrease with time.
(a)
(b)
(c)
Figure 13 shows the transient PDF of the leak-off velocity characterized by a normal distribution function. The leak-off velocity is the Darcy flow at the shale-fracture interface, which depends on the pore pressure gradient and is highest in the beginning of the injection (m s-1) and slows down to (m s-1) after 30 (min) of injection, as the pressure gradient nears zero.
Finally, Figure 14 shows the PDF of the transient fracture width, assuming the variability of the four input parameters (Table 4) for time interval [0,60] (min).
The mean fracture width is (mm) at and (mm) at (min). All possible values of fracture width fit a normal distribution and fall within a narrow range with a standard deviation of (mm) at , then the range of variation widens, and the standard deviation is (mm) at (min). The reason for this is that, at early times, the fracture width evolution is an elastic problem, depending primarily on one parameter Young’s modulus , and with time, becomes a porothermoelastic problem that depends on both the thermal and hydraulic properties, therefore simultaneously considering random parameters (Young’s modulus , linear thermal expansion coefficient the mobility ratio , and the thermal conductivity ).
4. Discussion of Findings
When ranking the parameters by order of influence on the solution, one must pay attention to the parameters with high variability as this may affect the sensitivity score. In this problem, the mobility ratio, thermal conductivity, and thermal expansion coefficient are high variability parameters (Table 3) and are also the dominating parameters for the pressure solution (Figure 7). Although high variability affects the sensitivity score, it does not eliminate the reliability of the method, because its main purpose is to rank the contributing parameters, as can be seen in Figure 8 where does not have a large variability but is the second most important parameter in the temperature solution. Results also show that temperature takes longer than pressure to stabilize in the system, clearly depicted in the resonance of the sensitivity curves in Figure 8. While at 30 (min), the pressure sensitivity ranking curves stabilize, and the temperature ones stabilize only at 60 (min). Figure 9 shows that the sensitivity ranking for the leak-off velocity, that is, Darcy flow at the shale-fracture interface, depends primarily on the mobility ratio, where its sensitivity score remains constant after 5 (min), thus indicating that the pressure gradient tends to zero and consequently generates a slow leak-off. Figure 10 shows that the fracture width evolution depends initially only on the elastic parameter and, with time, becomes mainly driven by thermal effects and depends primarily on the thermal expansion coefficient of the rock. Figure 11 shows that the applied excess pore pressure at the shale-fracture interface leads to an instantaneous pore pressure jump in the medium with a variability of [30-35] (MPa), and then, as time progresses, pressure diffuses further and deeper into the formation. We can also note that the normal distribution of the PDF of the pressure becomes narrower with time, meaning that as time progresses, the effects of the variability in formation properties decrease. Figure 12 clearly shows the temperature diffusion, where the cooling front reaches (m) after 15 (min) and (m) at 30 (min) and (m) at 60 (min). Results also show that the standard deviation of the temperature PDF is independent of time. Further, Darcy flow at the shale-fracture interface (the leak-off) stabilizes after 20 (min) (Figure 13), thus meaning that pressure at the boundary is reaching a constant value. Finally, the PDF of the fracture width provides its range of possible values and corresponding likelihood, which in turn determines if proppants will fit inside the created fracture. Figure 15 shows the “most likely” fracture width for three different time steps, Because of the thermal gradient at the shale-fracture interface, fracture width is continuously widening to reach a range of 0-1.5 (mm) after 60 (min) of injection, which means that one must wait before pumping proppants to profit from this mechanism. [3] found that fracture aperture can be as much as 70% larger when the temperature gradient at the interface is 60 (°C), contrary to when the rock and the injected fluid have the same temperature, in which case the fracture width starts to decrease shortly after pumping.
A proper combination of fracturing fluid and proppant kind/mesh size is paramount for a successful fracture job. The industry designs proppants to migrate far enough, thus maintaining a sufficient length/width of the fracture, and to keep the induced hydraulic fracture open, thus providing good fracture conductivity. According to API standards, two types of proppant mesh/size are primarily used in the field: with a diameter of 0.69 (mm) and with a diameter of 0.33 (mm). I test the agreement between the width PDF and the industry preferred proppant mixture of ( (mm)) (60%) and (40%) and find that the probability that the proppants and will fit inside the fracture at is 63% and 100%, respectively (Figure 16). Consequently, the mixture of 60% and 40% has a probability success rate of 100%, which is observed in the field.
Figure 17 shows the PDF for the fracture aperture at (min), where the probability that the proppants and will fit the fracture is 68% and 96%, respectively. This explains the success rate of the 60% and 40% mixture even at larger times.
5. Conclusions
I investigate here two important parts of uncertainty quantification in hydraulic fracture responses: sensitivity analysis and probability density functions. The study examines the impact of the uncertainty in the formation properties on the one-dimensional problem output solution. The analysis indicates that a large variability in the formation properties propagates a large uncertainty in the fracture responses, most importantly the fracture width. The result is consistent with field observations of successful proppant size, thus proving that uncertainty quantification for fracture width can be a powerful decision-making tool to seek critical proppant mesh/size.
Nevertheless, the investigated problem has some limitations. First, the one-dimensional solution is limited to model the flow of a Newtonian fluid with a constant viscosity, independent from both pressure and temperature. Further, we assume the fracturing fluid pressure and temperature inside the apertures to be constant throughout the stimulation period and only address the leak-off and heat transfer from the fracture into the formation.
Secondly, despite the simplicity of the Monte-Carlo methods and their robustness, they have the disadvantage of a slow convergence rate and consequently require a large sampling of the stochastic (or random) input parameters. To determine the variability induced on the solution, one must repeatedly calculate the solution with the generated population of random input parameters. Furthermore, these methods are limited to the local variability of the solution because the impact of the input uncertainties is not simultaneously considered (one at a time) to avoid long computation times and computer crashes.
In conclusion, uncertainty quantification (UQ) is a large field of study where some methods may be more effective than others depending on the complexity of the problem. The choice of the optimum statistical method is usually achieved through trial and error and past literature review. Other statistical methods characterized by lower computational time and population size include the spectral methods and dimensional analysis and are worth further developing with application to hydraulic fracturing.
Data Availability
Data is available upon reasonable request.
Conflicts of Interest
The author declares that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The author appreciates helpful comments from Stephen A. Miller. This publication is supported by the Université de Neuchâtel.