# Uncertainty Quantification for a Hydraulic Fracture Geometry: Application to Woodford Shale Data

Authors: Batoul M. Gisler
Journal: Geofluids
 (2021)

---

## 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.


---

## Body

## 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).Figure 1
Simplified drawing of the natural fractures in Woodford shale.Figure 2
Conceptual model for hydraulic fracturing in shale reservoir with preexisting natural fractures. Horizontal well injection initiates multibranch fractures. I delimit the study to one multibranch fracture highlighted with a purple dashed box.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.Figure3 zooms in to a multibranch fracture and shows how the fracture apertures are perpendicular to the minimum horizontal stress σh,min and intersect a porous medium characterized with an initial pore pressure p0 and an original temperature T0 and subject to fracturing fluid with pressure pf and temperature Tf at its faces (Figure 3). Figure 3 also shows that the distance between two consecutive fracture apertures is 2h, but because of the mirror symmetry, the one-dimensional solution is constrained to the length h. 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.Figure 3
The red dashed box zooms in to two consecutive fractures and shows the fracture geometry with the initial and boundary conditions:Pf and Tf are the fluid pressure and temperature, respectively. P0 and T0 are the initial pressure and temperature, respectively. The 1D solution is limited to h because of the 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 (xy) plane is the isotropic plane and the z 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.Figure 4
One-dimensional consolidation setup showing the initial and top boundary conditions, where the plane (xy) is the isotropic plane, where the x axis has a subscript of 1, and the z axis is the axis of rotational symmetry with a subscript of 3. The top, bottom, and right boundaries (in bold) are assumed impermeable, adiabatic, and frictionless.First, the analytical pressure solution is derived as(1)px,t=2h∑m=0∞1ξc1e−γ1ξ2t+c2e−γ2ξ2tsinξx+ΔP,then, the temperature is given by
(2)Tx,t=2h∑m=0∞1ξc3e−γ1ξ2t+c4e−γ2ξ2tsinξx+ΔT,finally, the displacement is given by
(3)ux,t=S0−α1ΔP−β1sΔTM11x−h+2α1M11h∑m=0∞1ξ2c1e−γ1ξ2t+c2e−γ2ξ2tcosξx+2β1sM11h∑m=0∞1ξ2c3e−γ1ξ2t+c4e−γ2ξ2tcosξx,where h is the half distance (m), ξ is 2m+1π/2h (m-1), S0 is the stress boundary condition that is the difference between the fracturing fluid pressure Pf (MPa) and the minimum horizontal stress Shmin (MPa), M11 is the stiffness coefficient (MPa) calculated from Young’s modulus and Poisson’s ratio, β1s is a thermal parameter calculated from the stiffness coefficients and the solid thermal expansion coefficient (MPa°C−1) and α1 is Biot’s coefficient [1], and c1,c2,c3,c4, γ1, and γ2 are lumped parameters defined as
(4)c1=γ1γ2g11p0−ΔP−g12ΔT−p0−ΔPγ1−γ1+γ2;c2=p0−ΔP−c1;c3=c11−γ1g11γ1g12;c4=−ΔT−c3;γ1=g11+g22+g11+g222−4g∗2g∗;γ2=g11+g22−g11+g222−4g∗2g∗;and g11,g12,g21,g22, and g∗ are defined as
(5)g11=μk11M+α12M11;g12=μk1α1β1sM11−βsf;g21=α1β1sT0λ1M11;g22=1λ1ρCv+β1s2T0M11;g∗=g22g11−g12g21,where μ is the fluid viscosity (Pa s), k1 is the permeability (m2), M is a poroelastic parameter in (MPa), λ1 is the thermal conductivity Jm−1s−1°C−1, and ρCv is the overall thermal capacity (J/m3°C). Note that subscripts 1 and 3 represent the x direction and the z direction, respectively.The initial and boundary conditions associated with Equations (1)–(3) are given as follows:
(6)px,0=p0;p0,t=pf;Tx,0=T0;T0,t=Tf;∂p∂xh,t=0;∂T∂xh,t=0;uh,t=0.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 byE, 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], M Biot’s modulus (GPa), and κ the mobility ratio (ratio between permeability of the formation and viscosity of the pore fluid: k/μ) (m2Pa−1s−1). 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 αis°C−1 in the ith direction, the volumetric thermal expansion coefficient of the pore fluid αf°C−1, the thermal conductivity of the porous medium λ1 (in the x direction) Jm−1s−1°C−1, and the overall thermal capacity ρCvJm−3°C−1. 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.Table 1
Formation properties used to model the uncertainty in the hydraulic fracture responses.
EiYoung’s modulus in theith directionνiPoisson’s ratioαisSolid skeleton linear thermal expansion coefficient in theith directionκ1Mobility ratio in thex directionλ1Thermal conductivity in thex directionρCvOverall thermal capacityIn 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 d1,d2,d3,⋯, each characterized by a known probability distribution law, and calculates the correspondent solution s of the model s1,s2,s3,⋯, 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 M−1/2 governs all Monte-Carlo methods, thus requiring a large population size M [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,y=fxi=1,xi=2,⋯.xi=m with m input parameters.
(1)
Considerx1 a random input parameter that fits a defined distribution function, e.g., normal distribution, with an assigned mean and standard deviation, and generate Npopulation of x1. Figure 5 shows the population distribution for Young’s modulus E1 (GPa) converging for Npopulation=15,000(2)
Calculate the deterministic solution for each value ofx1(3)
Calculate the variance of theNpopulation values of output solution and note it as y1(4)
Repeat procedure for each of the other input parametersi=2:m(5)
Calculate the sum of the output solutionyi for m values of variance. ∑i=0myi(6)
Calculate the coefficient of variationy1/∑i=1myi for each input parameter. This coefficient represents the degree of influence of each input parameter on the output solutionFigure 5
Monte-Carlo population simulation for Young’s modulus withNpopulation=15,000 with an average of 8.21 (GPa) and a standard deviation of 1.2 (GPa) assuming a normal distribution.


### 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 solutiony=fxi=1,xi=2,⋯.xi=m with m 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 ofn values for each of the selected random variables and calculate the output solution(3)
Plot the histogram of all the n output valuesy(4)
Fit the histogram to a distribution function and calculate the PDF function


## 2.1. Conceptual Model

Analytical solutions for pressure, temperature, strain, and displacement are derived in detail by [3] following the 1D consolidation theory.Figure3 zooms in to a multibranch fracture and shows how the fracture apertures are perpendicular to the minimum horizontal stress σh,min and intersect a porous medium characterized with an initial pore pressure p0 and an original temperature T0 and subject to fracturing fluid with pressure pf and temperature Tf at its faces (Figure 3). Figure 3 also shows that the distance between two consecutive fracture apertures is 2h, but because of the mirror symmetry, the one-dimensional solution is constrained to the length h. 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.Figure 3
The red dashed box zooms in to two consecutive fractures and shows the fracture geometry with the initial and boundary conditions:Pf and Tf are the fluid pressure and temperature, respectively. P0 and T0 are the initial pressure and temperature, respectively. The 1D solution is limited to h because of the 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 (xy) plane is the isotropic plane and the z 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.Figure 4
One-dimensional consolidation setup showing the initial and top boundary conditions, where the plane (xy) is the isotropic plane, where the x axis has a subscript of 1, and the z axis is the axis of rotational symmetry with a subscript of 3. The top, bottom, and right boundaries (in bold) are assumed impermeable, adiabatic, and frictionless.First, the analytical pressure solution is derived as(1)px,t=2h∑m=0∞1ξc1e−γ1ξ2t+c2e−γ2ξ2tsinξx+ΔP,then, the temperature is given by
(2)Tx,t=2h∑m=0∞1ξc3e−γ1ξ2t+c4e−γ2ξ2tsinξx+ΔT,finally, the displacement is given by
(3)ux,t=S0−α1ΔP−β1sΔTM11x−h+2α1M11h∑m=0∞1ξ2c1e−γ1ξ2t+c2e−γ2ξ2tcosξx+2β1sM11h∑m=0∞1ξ2c3e−γ1ξ2t+c4e−γ2ξ2tcosξx,where h is the half distance (m), ξ is 2m+1π/2h (m-1), S0 is the stress boundary condition that is the difference between the fracturing fluid pressure Pf (MPa) and the minimum horizontal stress Shmin (MPa), M11 is the stiffness coefficient (MPa) calculated from Young’s modulus and Poisson’s ratio, β1s is a thermal parameter calculated from the stiffness coefficients and the solid thermal expansion coefficient (MPa°C−1) and α1 is Biot’s coefficient [1], and c1,c2,c3,c4, γ1, and γ2 are lumped parameters defined as
(4)c1=γ1γ2g11p0−ΔP−g12ΔT−p0−ΔPγ1−γ1+γ2;c2=p0−ΔP−c1;c3=c11−γ1g11γ1g12;c4=−ΔT−c3;γ1=g11+g22+g11+g222−4g∗2g∗;γ2=g11+g22−g11+g222−4g∗2g∗;and g11,g12,g21,g22, and g∗ are defined as
(5)g11=μk11M+α12M11;g12=μk1α1β1sM11−βsf;g21=α1β1sT0λ1M11;g22=1λ1ρCv+β1s2T0M11;g∗=g22g11−g12g21,where μ is the fluid viscosity (Pa s), k1 is the permeability (m2), M is a poroelastic parameter in (MPa), λ1 is the thermal conductivity Jm−1s−1°C−1, and ρCv is the overall thermal capacity (J/m3°C). Note that subscripts 1 and 3 represent the x direction and the z direction, respectively.The initial and boundary conditions associated with Equations (1)–(3) are given as follows:
(6)px,0=p0;p0,t=pf;Tx,0=T0;T0,t=Tf;∂p∂xh,t=0;∂T∂xh,t=0;uh,t=0.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 byE, 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], M Biot’s modulus (GPa), and κ the mobility ratio (ratio between permeability of the formation and viscosity of the pore fluid: k/μ) (m2Pa−1s−1). 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 αis°C−1 in the ith direction, the volumetric thermal expansion coefficient of the pore fluid αf°C−1, the thermal conductivity of the porous medium λ1 (in the x direction) Jm−1s−1°C−1, and the overall thermal capacity ρCvJm−3°C−1. 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.Table 1
Formation properties used to model the uncertainty in the hydraulic fracture responses.
EiYoung’s modulus in theith directionνiPoisson’s ratioαisSolid skeleton linear thermal expansion coefficient in theith directionκ1Mobility ratio in thex directionλ1Thermal conductivity in thex directionρCvOverall thermal capacityIn 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 d1,d2,d3,⋯, each characterized by a known probability distribution law, and calculates the correspondent solution s of the model s1,s2,s3,⋯, 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 M−1/2 governs all Monte-Carlo methods, thus requiring a large population size M [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,y=fxi=1,xi=2,⋯.xi=m with m input parameters.
(1)
Considerx1 a random input parameter that fits a defined distribution function, e.g., normal distribution, with an assigned mean and standard deviation, and generate Npopulation of x1. Figure 5 shows the population distribution for Young’s modulus E1 (GPa) converging for Npopulation=15,000(2)
Calculate the deterministic solution for each value ofx1(3)
Calculate the variance of theNpopulation values of output solution and note it as y1(4)
Repeat procedure for each of the other input parametersi=2:m(5)
Calculate the sum of the output solutionyi for m values of variance. ∑i=0myi(6)
Calculate the coefficient of variationy1/∑i=1myi for each input parameter. This coefficient represents the degree of influence of each input parameter on the output solutionFigure 5
Monte-Carlo population simulation for Young’s modulus withNpopulation=15,000 with an average of 8.21 (GPa) and a standard deviation of 1.2 (GPa) assuming a normal distribution.


## 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 solutiony=fxi=1,xi=2,⋯.xi=m with m 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 ofn values for each of the selected random variables and calculate the output solution(3)
Plot the histogram of all the n output valuesy(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 Table2 [17–20]. The transversely isotropic nature of this low permeability shale signifies that its properties depend on the direction, where x has a subscript of 1 and z 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 αs. The model considers thermal effects because of the significant temperature difference ∆T=60°C between the fracturing fluid temperature Tf=30°C and the formation temperatureT0=90°C.Table 2
Woodford shale data used to test the sensitivity of the model to the input uncertainties.
Mobility ratioκ1,m2/Pas45Solid skeleton linear thermal expansion coefficientαs, /°C2×10−5Overall heat capacityρCv, J/m3°C2.75Thermal conductivityλ1,J/ms°C1.4Original pore pressurepo, MPa27.6Original temperatureTo, °C90Minimum horizontal stressσh,min, MPa52Stress boundary conditionS0, MPa12Fracturing fluid pressurepf, MPa64Fracturing fluid temperatureTf, °C30Young’s modulusE1,E3, GPa5.6, 4.6Poisson’s ratiosν1,ν30.3Distance between vertical apertures2h, m1.2However, literature review [17, 19, 20] confirms a large variability in some of those properties, summarized in Table 3.Table 3
Variability of some Woodford shale formation properties from field measurements.
Mobility ratioκ1, nD/cP20-60Solid skeleton linear thermal expansion coefficientαs, /°C0.9−4.8×10−5Overall thermal capacityρCv, J/m3°C2.57-2.94Thermal conductivityλ1, J/ms°C0.7-1.4Young’s modulusE1,E3,GPa6.8-9.8; 3.85-5.6Poisson’s ratioν1,ν30.1-0.15; 0.23-0.3


### 3.2. Sensitivity Analysis: Application to Woodford Shale

The saturated porous medium in Figure4 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 Pf, and consequently, the sensitivity score for all parameters is zero at x=0 (m) (Figure 6). Figure 6 shows that the dominant control on the initial pore pressure (Equation (1)) is the elastic Young’s modulus E1 with a sensitivity score of 90%. As pressure distribution stabilizes inside the half-space, the sensitivity of the pore pressure to E1 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 E1 significantly decreases and the pore pressure becomes more sensitive to the mobility ratio κ1, the thermal conductivity λ1, and the skeleton linear thermal expansion coefficient αs, 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.Figure 6
Sensitivity ranking for pressure in the fracture at timet=0+. The instantaneous pore pressure jump is 90% governed by Young’s modulus E1.Figure 7
Sensitivity ranking for pressure in the fracture at time (a)t=15 (min), (b) t=30 (min), and (c) t=60 (min).
(a)(b)(c)λ1 is the second most dominant parameter, where its influence reaches 0.09 (m) at t=15 (min) (Figure 7(a)), 0.12 (m) at t=30 (min) (Figure 7(b)), and 0.18 (m) at t=60 (min) (Figure 7(c)), and αs is the third most dominant parameter, where its influence reaches 0.06 (m) at t=15 (min) (Figure 7(a)), 0.09 (m) at t=30 (min) (Figure 7(b)), and 0.12 (m) at t=60 (min) (Figure 7(c)). The thermal conductivity λ1 sensitivity curve presents a peak that coincides with the dip in the mobility ratio sensitivity curve for x=0.04 (m) at t=15 (min) (Figure 7(a)), x=0.05 (m) at t=30 (min) (Figure 7(b)), and x=0.06 [m] at t=60 (min) (Figure 7(c)), highlighting a lower pressure zone, as seen in [3]. After 30 (min) of injection, the dip of the mobility ratio κ1 stabilizes at 48% and the peak in the thermal conductivity λ1 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 E1 and κ1, such as x>0.1 (m) at t=15 (min) (Figure 7(a)), x>0.12 (m) at t=30 (min) (Figure 7(b)), and x>0.18 (m at t=60 (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 t=0+, because the temperature front takes finite time to diffuse inside the formation.Figure8(a) shows that at t=15 (min), the cooling front reaches 0.11 (m) and depends on the thermal conductivity λ1 and the overall thermal capacity ρCv. Figures 8(b) and 8(c) show that the sensitivity of the temperature to ρCv increases while the influence of λ1 continues to decrease but remains dominant, with a sensitivity score of 70% for λ1 and for of 30% for ρCv at t=60 (min) and for x=0.2 (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 (x=0 (m)) only depends on the thermal conductivity λ1. Note that the Monte-Carlo simulation for temperature becomes smoother as the system stabilizes and reaches a steady condition.Figure 8
Sensitivity ranking for temperature in the fracture at time (a)t=15 (min), (b) t=30 (min), and (c) t=60 (min).
(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 (x=0 (m)).Figure9 shows that the leak-off velocity is most sensitive to the mobility ratio κ1 and the thermal conductivity λ1, with a sensitivity score of 90% and 10%, respectively.Figure 9
Sensitivity ranking for the transient leak-off velocity (m/s) for time intervalt=0,60 (min).Finally, the fracture width is the displacement at the shale-fracture interface,x=0 (m) (Equation (3)).Similar to the pressure (Figure6(a)), fracture width is initially (t=0+) most sensitive to the elastic Young’s modulus E1 (Figure 10), but because thermal effects are factored in the rock deformation process, the influence of E1 decreases with time, while the linear thermal expansion coefficient αs, the mobility ratio κ1, and the thermal conductivity λ1 start to play a larger role in the displacement solution. Figure 10 shows that after 60 (min), the fracture width becomes most sensitive to αs (75%), followed by κ1 (12%), λ1 (9%), and E1 (3%).Figure 10
Sensitivity ranking for the fracture width (mm) for time intervalt=0,60 (min).Table4 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.Table 4
Sensitivity analysis summary showing which input parameters influence the solution. Y is for yes and N is for No.
OutputParametersαsE1E3ν1ν3λ1ρCvκ1PressureYYNNNYNYTemperatureNNNNNYYNLeak-offNNNNNNNYFracture widthYYNNNYNY


### 3.3. Probability Density Functions: Application to Woodford Shale

Figures11(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 t=15 (min), t=30 (min), and t=60 (min), respectively. The pore pressure variability fits a normal distribution function, where the pressure is fixed at 64 (MPa) at the shale-fracture interface, x=0 (m). That is, PDF at this boundary is 100% and is not included in the plots.Figure 11
Probability density functions for pressure at times (a)t=15 (min), (b) t=30 (min), and (c) t=60 (min).
(a)(b)(c)Figure11(a) shows that the transient pore pressure front diffuses from the boundary and reaches a depth of x=0.1 (m) at t=15 (min), to stabilize at a pressure range of 30-35 (MPa), for x>0.1 (m), which corresponds to the pore pressure jump felt in the porous medium at t=0+. 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 t=15 (min), which decreases with time, where the PDF curve becomes sharper at 30 (min) and 60 (min).Figures12(a)–12(c) show the probability density functions (PDF) of the thermal diffusion, assuming a random thermal conductivity λ1 and a random overall heat capacity ρCv, calculated at times t=15 (min), t=30 (min), and t=60 (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 t=60 (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.Figure 12
Probability density functions for temperature at times (a)t=15 (min), (b) t=30 (min), and (c) t=60 (min).
(a)(b)(c)Figure13 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 1.5×10−3 (m s-1) and slows down to 0.3×10−3 (m s-1) after 30 (min) of injection, as the pressure gradient nears zero.Figure 13
Probability density function of the leak-off velocity for time intervalt=0,60 (min).Finally, Figure14 shows the PDF of the transient fracture width, assuming the variability of the four input parameters (Table 4) for time interval [0,60] (min).Figure 14
Probability density function for fracture width for time intervalt=0,60 (min).The mean fracture width is0.72 (mm) at t=0+ and 0.84 (mm) at t=60 (min). All possible values of fracture width fit a normal distribution and fall within a narrow range with a standard deviation of σ=0.04 (mm) at t=0+, then the range of variation widens, and the standard deviation is σ=0.15 (mm) at t=60 (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 E1, and with time, becomes a porothermoelastic problem that depends on both the thermal and hydraulic properties, therefore simultaneously considering random parameters (Young’s modulus E1, linear thermal expansion coefficient αs, the mobility ratio κ1, and the thermal conductivity λ1).


## 3.1. Woodford Shale Data

I test the proposed framework on the 1D solution with typical Woodford Shale data presented in Table2 [17–20]. The transversely isotropic nature of this low permeability shale signifies that its properties depend on the direction, where x has a subscript of 1 and z 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 αs. The model considers thermal effects because of the significant temperature difference ∆T=60°C between the fracturing fluid temperature Tf=30°C and the formation temperatureT0=90°C.Table 2
Woodford shale data used to test the sensitivity of the model to the input uncertainties.
Mobility ratioκ1,m2/Pas45Solid skeleton linear thermal expansion coefficientαs, /°C2×10−5Overall heat capacityρCv, J/m3°C2.75Thermal conductivityλ1,J/ms°C1.4Original pore pressurepo, MPa27.6Original temperatureTo, °C90Minimum horizontal stressσh,min, MPa52Stress boundary conditionS0, MPa12Fracturing fluid pressurepf, MPa64Fracturing fluid temperatureTf, °C30Young’s modulusE1,E3, GPa5.6, 4.6Poisson’s ratiosν1,ν30.3Distance between vertical apertures2h, m1.2However, literature review [17, 19, 20] confirms a large variability in some of those properties, summarized in Table 3.Table 3
Variability of some Woodford shale formation properties from field measurements.
Mobility ratioκ1, nD/cP20-60Solid skeleton linear thermal expansion coefficientαs, /°C0.9−4.8×10−5Overall thermal capacityρCv, J/m3°C2.57-2.94Thermal conductivityλ1, J/ms°C0.7-1.4Young’s modulusE1,E3,GPa6.8-9.8; 3.85-5.6Poisson’s ratioν1,ν30.1-0.15; 0.23-0.3


## 3.2. Sensitivity Analysis: Application to Woodford Shale

The saturated porous medium in Figure4 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 Pf, and consequently, the sensitivity score for all parameters is zero at x=0 (m) (Figure 6). Figure 6 shows that the dominant control on the initial pore pressure (Equation (1)) is the elastic Young’s modulus E1 with a sensitivity score of 90%. As pressure distribution stabilizes inside the half-space, the sensitivity of the pore pressure to E1 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 E1 significantly decreases and the pore pressure becomes more sensitive to the mobility ratio κ1, the thermal conductivity λ1, and the skeleton linear thermal expansion coefficient αs, 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.Figure 6
Sensitivity ranking for pressure in the fracture at timet=0+. The instantaneous pore pressure jump is 90% governed by Young’s modulus E1.Figure 7
Sensitivity ranking for pressure in the fracture at time (a)t=15 (min), (b) t=30 (min), and (c) t=60 (min).
(a)(b)(c)λ1 is the second most dominant parameter, where its influence reaches 0.09 (m) at t=15 (min) (Figure 7(a)), 0.12 (m) at t=30 (min) (Figure 7(b)), and 0.18 (m) at t=60 (min) (Figure 7(c)), and αs is the third most dominant parameter, where its influence reaches 0.06 (m) at t=15 (min) (Figure 7(a)), 0.09 (m) at t=30 (min) (Figure 7(b)), and 0.12 (m) at t=60 (min) (Figure 7(c)). The thermal conductivity λ1 sensitivity curve presents a peak that coincides with the dip in the mobility ratio sensitivity curve for x=0.04 (m) at t=15 (min) (Figure 7(a)), x=0.05 (m) at t=30 (min) (Figure 7(b)), and x=0.06 [m] at t=60 (min) (Figure 7(c)), highlighting a lower pressure zone, as seen in [3]. After 30 (min) of injection, the dip of the mobility ratio κ1 stabilizes at 48% and the peak in the thermal conductivity λ1 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 E1 and κ1, such as x>0.1 (m) at t=15 (min) (Figure 7(a)), x>0.12 (m) at t=30 (min) (Figure 7(b)), and x>0.18 (m at t=60 (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 t=0+, because the temperature front takes finite time to diffuse inside the formation.Figure8(a) shows that at t=15 (min), the cooling front reaches 0.11 (m) and depends on the thermal conductivity λ1 and the overall thermal capacity ρCv. Figures 8(b) and 8(c) show that the sensitivity of the temperature to ρCv increases while the influence of λ1 continues to decrease but remains dominant, with a sensitivity score of 70% for λ1 and for of 30% for ρCv at t=60 (min) and for x=0.2 (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 (x=0 (m)) only depends on the thermal conductivity λ1. Note that the Monte-Carlo simulation for temperature becomes smoother as the system stabilizes and reaches a steady condition.Figure 8
Sensitivity ranking for temperature in the fracture at time (a)t=15 (min), (b) t=30 (min), and (c) t=60 (min).
(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 (x=0 (m)).Figure9 shows that the leak-off velocity is most sensitive to the mobility ratio κ1 and the thermal conductivity λ1, with a sensitivity score of 90% and 10%, respectively.Figure 9
Sensitivity ranking for the transient leak-off velocity (m/s) for time intervalt=0,60 (min).Finally, the fracture width is the displacement at the shale-fracture interface,x=0 (m) (Equation (3)).Similar to the pressure (Figure6(a)), fracture width is initially (t=0+) most sensitive to the elastic Young’s modulus E1 (Figure 10), but because thermal effects are factored in the rock deformation process, the influence of E1 decreases with time, while the linear thermal expansion coefficient αs, the mobility ratio κ1, and the thermal conductivity λ1 start to play a larger role in the displacement solution. Figure 10 shows that after 60 (min), the fracture width becomes most sensitive to αs (75%), followed by κ1 (12%), λ1 (9%), and E1 (3%).Figure 10
Sensitivity ranking for the fracture width (mm) for time intervalt=0,60 (min).Table4 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.Table 4
Sensitivity analysis summary showing which input parameters influence the solution. Y is for yes and N is for No.
OutputParametersαsE1E3ν1ν3λ1ρCvκ1PressureYYNNNYNYTemperatureNNNNNYYNLeak-offNNNNNNNYFracture widthYYNNNYNY


## 3.3. Probability Density Functions: Application to Woodford Shale

Figures11(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 t=15 (min), t=30 (min), and t=60 (min), respectively. The pore pressure variability fits a normal distribution function, where the pressure is fixed at 64 (MPa) at the shale-fracture interface, x=0 (m). That is, PDF at this boundary is 100% and is not included in the plots.Figure 11
Probability density functions for pressure at times (a)t=15 (min), (b) t=30 (min), and (c) t=60 (min).
(a)(b)(c)Figure11(a) shows that the transient pore pressure front diffuses from the boundary and reaches a depth of x=0.1 (m) at t=15 (min), to stabilize at a pressure range of 30-35 (MPa), for x>0.1 (m), which corresponds to the pore pressure jump felt in the porous medium at t=0+. 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 t=15 (min), which decreases with time, where the PDF curve becomes sharper at 30 (min) and 60 (min).Figures12(a)–12(c) show the probability density functions (PDF) of the thermal diffusion, assuming a random thermal conductivity λ1 and a random overall heat capacity ρCv, calculated at times t=15 (min), t=30 (min), and t=60 (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 t=60 (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.Figure 12
Probability density functions for temperature at times (a)t=15 (min), (b) t=30 (min), and (c) t=60 (min).
(a)(b)(c)Figure13 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 1.5×10−3 (m s-1) and slows down to 0.3×10−3 (m s-1) after 30 (min) of injection, as the pressure gradient nears zero.Figure 13
Probability density function of the leak-off velocity for time intervalt=0,60 (min).Finally, Figure14 shows the PDF of the transient fracture width, assuming the variability of the four input parameters (Table 4) for time interval [0,60] (min).Figure 14
Probability density function for fracture width for time intervalt=0,60 (min).The mean fracture width is0.72 (mm) at t=0+ and 0.84 (mm) at t=60 (min). All possible values of fracture width fit a normal distribution and fall within a narrow range with a standard deviation of σ=0.04 (mm) at t=0+, then the range of variation widens, and the standard deviation is σ=0.15 (mm) at t=60 (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 E1, and with time, becomes a porothermoelastic problem that depends on both the thermal and hydraulic properties, therefore simultaneously considering random parameters (Young’s modulus E1, linear thermal expansion coefficient αs, the mobility ratio κ1, and the thermal conductivity λ1).


## 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 (Table3) 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 ρCv 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 E1 and, with time, becomes mainly driven by thermal effects and depends primarily on the thermal expansion coefficient αs 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 x<0.08 (m) after 15 (min) and x<0.12 (m) at 30 (min) and x<0.175 (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, t=0+,t=30min,andt=60min. 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.Figure 15
Probability density function of the fracture width at different time steps.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:D20/40 with a diameter of 0.69 (mm) and D40/70 with a diameter of 0.33 (mm). I test the agreement between the width PDF and the industry preferred proppant mixture of D100 (0.152 (mm)) (60%) and D40/70 (40%) and find that the probability that the proppants D20/40 and D40/70 will fit inside the fracture at t=0+ is 63% and 100%, respectively (Figure 16). Consequently, the mixture of 60% D100 and 40% D40/70 has a probability success rate of 100%, which is observed in the field.Figure 16
PDF of the fracture width att=0+. The probability of the proppant to fit inside the created fracture width is represented by the hashed area. The success rate of 63% for D20/40 and of 100% for D40/70.Figure17 shows the PDF for the fracture aperture at t=60 (min), where the probability that the proppants D20/40 and D40/70 will fit the fracture is 68% and 96%, respectively. This explains the success rate of the 60% D100 and 40% D40/70 mixture even at larger times.Figure 17
Probability density function of the fracture width att=60 (min). The probability of the proppant to fit inside the created fracture width is represented by the hashed area. The success rate of 68% for D20/40 and 96% for D40/70.


## 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.



---

*Source: 2138115-2021-08-23.xml*