Sampling-based stochastic analysis of the PKN model for hydraulic fracturing

Hydraulic fracturing processes are surrounded by uncertainty, as available data is typically scant. In this work, we present a sampling-based stochastic analysis of the hydraulic fracturing process by considering various system parameters to be random. Our analysis is based on the Perkins-Kern-Nordgren (PKN) model for hydraulic fracturing. This baseline model enables computation of high fidelity solutions, which avoids pollution of our stochastic results by inaccuracies in the deterministic solution procedure. In order to obtain the desired degree of accuracy of the computed solution, we supplement the employed time-dependent moving-mesh finite element method with two new enhancements: (i) global conservation of volume is enforced through a Lagrange multiplier; (ii) the weakly singular behavior of the solution at the fracture tip is resolved by supplementing the solution space with a tip enrichment function. This tip enrichment function enables the computation of the tip speed directly from its associated solution coefficient. A novel incremental-iterative solution procedure based on a backward-Euler time-integrator with sub-iterations is employed to solve the PKN model. Direct Monte-Carlo sampling is performed based on random variable and random field input parameters. The presented stochastic results quantify the dependence of the fracture evolution process—in particular the fracture length and fracture opening—on variations in the elastic properties and leak-off coefficient of the formation, and the height of the fracture.


Introduction
Hydraulic fracturing processes -i.e. the fracturing of rock formations by a pressurized liquid to improve connectivity of reservoirs and geothermal formations -are surrounded by uncertainty, as available formation data is limited. Further improvement of models and simulation tools to understand this process is instrumental to increasing operational effectivity and to reliably quantify the risks and uncertainties that are involved.
Hydraulic fracturing models involve the coupling of three sub-models: i) a solid mechanics model which describes the deformation of the rock formation induced by the fluid load; ii) a fluid flow model to describe the fracturing-fluid in the crack, as well as its leak-off into the rock formation; iii) a fracture mechanics model including a fracture propagation criterion. Intrinsic characteristics of such models are the (strong) non-linearities related to the coupling between the solid and fluid, the singularities in the physical fields near the fracture front, the moving (fracture) domain boundaries, the degeneration of the governing equations near the tip region, and pronounced multiscale effects with the length scales varying from millimeters for the fracture opening near the tip to kilometers for the length of the fracture.
Various practical model simplifications -most commonly restricting the model to a single two-dimensional planar fracture -have been proposed, the most prominent of which are: i) The Perkins-Kern-Nordgren (PKN) model [1,2] for a fracture of fixed height and elliptical cross section (leveraging the Sneddon solution for the elasticity problem [3]) which propagates horizontally in one direction; ii) The radial model for a horizontal Penny-shaped crack that evolves evenly in all directions in accordance with Sneddon's solution [4]; iii) The Khristianovic-Geertsma-De Klerk (KGD) model [5,6,7] for a vertical fracture of fixed height that propagates horizontally in one direction. Various pseudo-three-dimensional (P3D) models [8] and planar-three-dimensional (PL3D) models [9,10] have been proposed over the years to enable consideration of more complex fracture patterns and fractures in multi-layer formations. These P3D models typically extend the above-mentioned two-dimensional models by considering a variation of fracture height in combination with fracture length and width. Although these classical models are generally based on restrictive and often ad-hoc assumptions, they are still widely used in the industry [11].
Despite the simplifications in the above-mentioned models, analytical solutions can only be obtained in limiting cases (see e.g. [5,6,2]). General solution strategies for these models rely on the use of computational techniques. Versatile and reliable simulators for hydraulic fracturing processes are indispensable in gaining further understanding of the process, in particular because direct observation possibilities (typically several kilometers below the earth surface) are limited. The most prominent computational techniques used for hydraulic fracturing are Finite Volume Methods [12], Finite Element Methods [13,14,15], Boundary Integral Methods [16] and Discrete Element Methods [17]. Recent advancements in numerical methods for hydraulic fracturing that are particularly noteworthy are the eXtended Finite Element Techniques (X-FEM) [18] and phasefield methods [19].
Although the importance of considering realistic geological situations is acknowledged [20], hydraulic fracture evolution in formations with uncertain heterogeneous rock properties (e.g. elastic moduli, compression/tensile strength, porosity, permeability) have not been studied in detail. The possibility of considering stochastic heterogeneities has been explored in [21], where the initiation and evolution of fracture has been studied using Monte Carlo sampling. The main challenge in such studies relates to the computational feasibility, in the sense that the computational effort of the deterministic simulations (whose error must be controlled in relation to the stochastic variations) makes them impractical in the context of sampling-based stochastic techniques.
In this work we present a detailed probabilistic analysis of the hydraulic fracturing process based on Monte Carlo simulations. The computational tractability of the stochastic framework considered herein motivates the use of the two-dimensional PKN model. We represent the (epistemic) uncertain parameters of the PKN model as random variables and/or random fields, and investigate the influence of these uncertain input parameters on the fracture geometry (in particular the fracture length and opening at the well bore). As part of the stochastic analysis we present a sensitivity analysis of the deterministic model. It is worth mentioning that similar analyses can i(t)  [22,23]. The primary focus of this work is the direct analysis of the propagation of heterogeneous uncertainties in the hydraulic fracturing process. As part of this study we present a detailed derivation of the PKN model in the context of random (spatially correlated) heterogeneous data. Control over accuracy is of paramount importance, as stochastic uncertainty quantification should not be polluted by numerical inaccuracies. In this regard we propose to use two features in the numerical method for the PKN model to control its accuracy: i) a Lagrange multiplier method to enforce the conservation of volume; ii) a special enrichment function for the finite element discretization of the PKN model to overcome tip singularity issues.
In Section 2 the governing equations are discussed, with a special focus on the incorporation of rock heterogeneities in the PKN model. In Section 3, the weak formulation and finite element discretization of the model are presented including various algorithmic details. The stochastic setting and Monte Carlo method are introduced in Section 4, where the random field discretization of the heterogeneous properties using the Karhunen-Loève expansion is also discussed. Numerical simulations are presented in Section 5 to study the influence of input uncertainties in hydraulic fracturing. Finally, conclusions are presented in Section 6.

The PKN model for hydraulic fracturing
In this section we review the PKN model for hydraulic fracturing in the context of stochastic analysis with heterogeneous random fields. The PKN model -which was originally formulated by Perkins and Kern [1] and later amended with a leak-off model by Nordgren [2] and a propagation condition by Kemp [24] -is a practical candidate for studying the probabilistic behavior of hydraulic fracturing by virtue of its computational tractability. Although highly simplified, the PKN model is based on fundamental physical principles and is capable of generating practically meaningful results [25].

Problem definition
The key assumption of the PKN model is that it considers a planar fracture with a constant height H; see Figure 1. Displacements and displacement gradients in the surrounding solid are assumed to remain small, and the material is assumed to be linear elastic and isotropic. The fracture surface resides in the xy-plane, while the fracture opens in the z-direction. The fracture aperture in the fracture plane is denoted by w(x, y, t), and the aperture at the y = 0 center line byŵ(x, t). The fracture connects to the well at x = 0 and its evolving front is situated at x = L(t).
A Newtonian fluid is injected into the fracture at the well with a controlled flow rate i(t), and the flow inside the permeable crack is assumed to be laminar. The fracture process is assumed to be in the viscosity-based regime, where toughness effects can be neglected (propagation is governed by friction and leak-off effects). At the front of the fracture fluid lag is assumed to be zero, i.e., the fracture front coincides with the fluid front. Moreover, spurt losses due to the creation of new fracture surface are ignored.

Governing equations
In this section the governing equations of the sub-models are reviewed. In the presented derivations we focus on those aspects of the sub-models that need careful consideration in the context of the stochastic analysis discussed in the remainder of this work.

Fluid flow model
The PKN model is based on the conservation of mass of the fluid, which establishes a link between the injected volume, the created fracture volume and the leak-off volume. The differential material balance for the fracturing fluid is given by where q(x, t) is the volume rate of flow through the cross-sectional area A(x, t) and s c l (x, t) is the rate of fluid volume loss per unit length of the fracture. At the well the flow rate is equal to q(0, t) = i(t).
The flow rate inside the fracture is related to the pressure gradient by assuming Poiseuille flow [26]. For such a flow the advective terms are assumed to remain small, so that the incompressible Navier-Stokes equations reduce to the Stokes equations. The PKN model moreover assumes a horizontal slot flow [27] with a parabolic fluid velocity profile where w(x, y, t) is the opening of the fracture and µ f is the fluid viscosity. As will be discussed in more detail in the context of the solid model, the assumptions of the PKN model render the cross-section to be ellipsoidal. The fracture aperture is then given by whereŵ is the maximum aperture at y = 0. The cross-sectional area is then equal to A(x, t) = π 4 Hŵ(x, t) and the fluid flow follows by integration of (2) as The leak-off volume rate per unit length of fracture in equation (1) follows the phenomenological law proposed by Carter [28]: In this expression c l is the leak-off coefficient and τ is the arrival time of the fracture tip at location L, i.e. τ = L −1 (x). The main assumptions behind this model are that: i) the fracturing fluid deposits a thin layer of relatively low permeability known as the filter cake on the inner faces of the fracture, with the deposition rate being proportional to the leak-off rate, and ii) the viscosity of the filtrate is high enough to fully displace the fluid already present in the rock pores. Substitution of equations (4) and (5) in the material balance (1) then yields the fluid flow mass balance:

Solid deformation model
To derive the relation between the fluid pressure and the solid deformation as used in the PKN model we consider the infinite domain Ω = R 3 + with evolving fracture surface Γ c (t) = {x ∈ Ω | 0 ≤ x ≤ L(t), −H/2 ≤ y ≤ y, z = 0}. Assuming inertia and gravity effects to be negligible, the solid deformation, u = (u, v, w), follows from the momentum balance in Ω where σ is the Cauchy stress. Along the fracture surface the solid is loaded by the fluid pressure, i.e. σn = −pn on Γ ± c where n is the normal vector pointing into the fluid domain. The Cauchy stress follows Hooke's law for isotropic materials where the Lamé parameters µ s (x) and λ s (x) are heterogeneous fields directly related to the Young's modulus, E(x), and Poisson's ratio, ν(x). In Section 4 we will model the elastic properties by means of random fields, which are parametrized by a mean value, a standard deviations, and an auto-correlation length, c. This auto-correlation length is a measure of the correlation between any two material points in a random field, where a large correlation length implies that the frequency of the heterogeneous field is low.
To derive the elasticity relation for the PKN model it is assumed that deformations are planar, in the sense that the solid does not deform in the direction of the fracture propagation (the x-direction), and a plane strain condition in that direction applies. The auto-correlation length must be sufficiently large not to violate these assumptions. Moreover, we herein neglect heterogeneous variations perpendicular to the fracture propagation direction. Under these assumptions, combination of the momentum balance (7) and constitutive relation (8) yields at an arbitrary plane perpendicular to the x-direction. Supplemented with the boundary conditions σ zz (x) = −p(x) and σ yz (x) = 0 on the fracture boundary and vanishing far field conditions, this problem can be solved analytically. The fracture opening in the case of a constant pressure in the yz-plane is given by (see e.g. Lowengrub [29] and Sneddon [3] for details) where is the plane strain modulus which is heterogeneous only in the x-direction. The elliptical profile (3) is a direct result of the setting of the elasticity problem considered here, with the maximum aperture equal tô An essential property of this solution is that along the crack path the fracture aperture is linearly related to the pressure. The local nature of this relation is a direct consequence of the assumed planar deformation. Note that the stress field and displacement field can be derived in the form of special integral representations [30].

Fracture propagation model
In the PKN model it is assumed that once the fracture has exceeded a certain distance, the energy dissipation associated with the fracture of the rock material is small compared to energy dissipation associated with the viscous flow of the fracturing fluid. This effectively neglects the fracture toughness, and fracture propagation is purely driven by the fluid velocity. Herein we adopt the standard assumption of zero fluid lag [2] -i.e., the velocity of the fluid at the fluid front and the tip propagation speed are equal -so that tip propagation follows the Stefan condition v(L(t), t) = lim Substitution of the flow rate (4) and surface area then yields:

The coupled initial boundary value problem
The hydraulic fracture problem is characterized by a strong coupling of the sub-models discussed above. The solid deformation is coupled to the fluid through the pressure loading along the fracture surface, while the fluid flow depends on the fracture opening though the Poiseuille flow profile. The fracture propagation condition is coupled directly to the fluid flow through the Stefan condition (12), and in turn influences the fluid flow by virtue of the fact that fracture propagation extends the fluid flow domain. The pointwise relation between pressure and the fracture opening in equation (11) allows the formulation of a single-field problem. Herein we consider the initial boundary value problem for the fracture opening on the time interval (0, T ) t with evolving domain (0, L(t)) x: ∀t ∈ (0, T ) Note that the omission of fluid lag in the model results in the tip boundary conditionŵ(L(t), t) = 0, reflecting zero fracture opening at the tip. This boundary condition leads to singular behavior of the fracture opening (and pressure) at the tip, which is an important characteristic of the coupled problem (14). The nature of this singularity is best observed from the tip-propagation relation (14e), where the plane strain modulus is assumed to be non-singular. Essentially, the propagation condition conveys that in order to have a finite propagation speed the derivative of the cube opening should be bounded. A function that reflects this behavior, while simultaneously satisfying the zero fracture opening tip condition, isŵ(x, t) ∝ 3 L(t) − x.

Deterministic computational methodology
In this section we present a computational methodology that enables the computation of solutions of the PKN model with an accuracy that makes it suitable for conducting a sampling-based stochastic analysis. In Section 3.1 we first discuss the incremental-iterative solution procedure which is used to integrate the time-dependent moving-boundary problem. In Section 3.2 we then discuss the spatial finite element discretization of the nonlinear system of equations introduced above, including two essential enhancements, viz. incorporation of a Lagrange multiplier to enforce the volume-conservation constraint, and a singular enrichment to resolve the singularity at the fracture tip.

Incremental-iterative solution procedure
To solve the time-dependent moving-boundary problem (14) we employ the incremental-iterative solution procedure outlined in Algorithm 1. We denote the time step size and index by ∆t and ı = 0, . . . , n t , respectively, such that t ı = ı∆t and T = n t ∆t. The solution at time step ı and sub-iteration  = 0, 1, 2, . . . is written as w ı  (x) and L ı  . The sub-iteration index  is omitted for converged solutions, i.e., w ı (x) and L ı .
We consider the implicit time-integration of both the fracture aperture and the fracture length, such that the coupled system (14) is discretized in time as Algorithm 1: Incremental-iterative solution procedure for all ı = 1, . . . , n t with initial conditionsŵ 0 (x) = 0 and L 0 = L 0 .
To solve this moving-boundary problem at time step ı = 1, . . . , n t , within each time step we sub-iterate between the aperture problem (15a)-(15c) and the propagation problem (15d) until convergence is achieved. The solution of the aperture problem (solve aperture in Algorithm 1) is approximated using a finite element discretization in combination with a Newton-Raphson procedure to resolve the non-linearity. Details regarding the finite element discretization will be discussed in Section 3.2. The propagation problem is solved by using the regula falsi method to find the root of the residual function Given an iterate for the fracture length, L ı  , the corresponding fracture aperture w ı  is computed using the solve aperture procedure, after which the procedure evaluate tip speed is called to determine the associated tip speed in accordance with: The regula falsi procedure is initialized with the fracture length at the previous time step, L ı 0 = L ı−1 , and with the forced propagation, L ı 1 = L ı−1 + ∆L. We note that the residual r ı 0 = −∆tL ı  is non-negative as a consequence of the non-negativity of the propagation speed. The residual r ı 1 = ∆L − ∆tL ı  is forced to be positive by selecting the numerical parameter ∆L > 0 to be sufficiently large. For all simulations in Section 5 we set ∆L equal to the element size at the fracture tip, and take a time-step that is sufficiently small to ensure positivity of the residual r ı 1 . The sub-iteration procedure is terminated when the fracture length converges to a specified tolerance, i.e., |L ı  − L ı −1 | < tol L . The converged iteratesŵ ı  (x) and L ı  for the fracture aperture and length are then used as the initial conditions for the next time step. The fracture arrival function τ ı (x) for the next time step is updated in the update tau routine. Since the arrival time function is evaluated by linear interpolation in the list {(Lî, tî)} î ı=0 , this update routine merely appends the converged fracture length to this list. Note that we assume zero leak-off in the initial crack, so that we do not need to evaluate the arrival function for x ≤ L 0 . For x > L ı−1 the arrival time is taken as the time at which the crack reached L ı−1 . As a result t ı − τ ı−1 (x) ≥ ∆t ∀x, which effectively resolves the occurrence of a singularity in the leak-off term at the fracture tip.

Finite element discretization
The solve aperture routine uses the finite element method to compute the fracture aperturê w ı  based on the fracture length L ı  and the aperture and arrival time at the previous time step, w ı−1 and τ ı−1 , respectively. To derive the finite element formulation the weak form of (15) is considered, where the sub-iteration index has been dropped for notational convenience: The time-dependent test and trial space V ı is defined such that the Dirichlet boundary conditions at the tip are satisfied, and all integrals in the above formulation are bounded. Note that the right-hand-side integral involving the fracture aperture at the previous time step is computed over the fracture domain at the previous time step, which reflects the fact that ahead of the fracture tip the aperture is equal to zero. Moreover, the initial crack is excluded from the integration domain for the right-hand-side term involving the leak-off, which results from the assumption that there is no leak-off in the initial crack. The weak formulation (18) is discretized using the finite element method by approximating the maximum aperture asŵ where the index set I ı contains the indices of the shape functions {N i } i∈I ı constructed over a mesh T ı that partitions the evolving domain (0, L ı ). The corresponding discrete trial and test space are given by The finite element discretization considered in this work is based on linear Lagrange finite elements, where the linear basis function associated with the tip node is constrained in accordance with the zero-aperture Dirichlet condition at the tip. Because of the nature of the solution we grade the mesh toward the tip by specifying the element size at the tip, the increase ratio between two adjacent elements, and the maximal element size that is approached toward the inflow boundary. A schematic representation of such a graded mesh is shown in Figure 3.
A discretization of (18) based on linear finite elements -even though graded toward the tip -leads to an unacceptable loss of accuracy at meshes that are computationally tractable within the scope of this work. This performance deterioration is a consequence of: i) the flux (4) being highly non-linearly dependent on the fracture aperture; and ii) the singular behavior at the tip not being represented by the linear finite element basis. In the remainder of this subsection we propose two method enhancements to ameliorate this performance degradation. In Section 5.1 the numerical performance of these enhancements will be assessed.

Mass conservation constraint
Although the weak formulation (18) is consistently derived from the mass balance equation (6), the finite element approximation of (18) violates both the local mass balance and its integrated global version. Since adequate representation of the conservation of mass is of critical importance for obtaining accurate solutions, we herein propose to enforce global conservation in our approximation. We obtain the global balance of mass by integration of the time-discrete version of (6) over the entire domain: This global balance clearly shows that the injected volume is conserved through leak-off through the fracture (first term), fracture opening (second term), and fracture propagation (third term). This global conservation law is enforced by multiplying it with a Lagrange multiplier, and appending it to the weak formulation (18).

Singular tip enrichment
As already briefly discussed in Section 2.3, the aperture solution to the problem (14) is singular at the tip as a consequence of the non-linear character of the coupled system of equations. In Ref. [31,32,11] it is shown that the aperture solution in the proximity of the tip is proportional to 3 L(t) − x, and that the tip velocity in accordance with equation (14e) is therefore finite. Evidently, this singular solution behavior is approximated poorly by the linear finite element basis. As a matter of fact, when solely using the linear finite element basis, the tip propagation speed will always vanish when evaluated through equation (14e). To improve the finite element approximation we enrich the test and trial space V ı with the tip asymptotic which we localize to the tip region using the partition-of-unity method [33]. The enriched finite element interpolation of the aperture is then given bŷ where the index set J ⊂ I contains the indices of the nodes that are enriched. In the numerical simulations considered in Section 5 we only enrich the linear finite element function associated with the tip, which we have found to be effective. A schematic representation of this enrichment is shown in Figure 3.

Stochastic setting
In this section we introduce the stochastic setting of the PKN problem. In Section 4.1 the Monte Carlo method which we employ in this work is briefly introduced, after which the parametrization of the random variables and random fields for the model parameters are discussed in Sec- In the remainder of this work we will consider some of the model parameters to be uncertain, viz. the plane strain modulus E , the leak-off coefficient, c l , and the fracture heightH. We use the tilde diacritic to indicate that these parameters are stochastic. As observable parameters we will focus on the final fracture length,x f (T ), and the maximum fracture mouth opening,w max (T ).

Direct Monte Carlo sampling
In this work we use direct Monte Carlo sampling to compute the stochastic response of the PKN model. The primary reason for using direct Monte Carlo sampling is that it does not pose any restrictions on the distributions of the model parameters and the observables. Moreover, the nonintrusive character of the method enables its direct application to the PKN model. More advanced stochastic techniques can aid in reducing the computational effort of the stochastic problem, but application of such techniques to the highly non-linear moving boundary problem considered here is non-standard and beyond the scope of the current work.
where the symbols µ d i and σ d i denote the mean and standard deviation, respectively. Evidently, the accuracy of the estimators depends on the number of samples N . Given a confidence level C µ for the estimated mean µ (omitting the subscript d i for notational convenience) -meaning that the estimated mean has a relative accuracy of ±(1 − C µ )/2 with probability C µ -the minimal number of samples can be estimated by [34,35] where Φ is the cumulative density function of a standard normal random variable, and V d i = σ d i /µ d i is the coefficient of variation of the random observabled i . A rough estimate of this coefficient of variation can be obtained using a Monte Carlo simulation with a small number of samples. From (24) it becomes apparent that a draw-back of the direct Monte Carlo method is the slow convergence of the sampling error (an increase in confidence level) with increase in the number of samples. In the context of computational models such as that considered in this work this practically means that the simulations are computationally very intensive in the case of practically meaningful confidence levels. The fact that Monte Carlo sampling is non-intrusive -in the sense that it is a method that does not interfere with the deterministic model -makes parallelization possible. We have implemented a parallel master/slave algorithm for our simulations, which shows excellent scalability.

Random variable and random field parametrization
In this work we represent the considered scalar model parameters,m i , by means of log-normal distributions, which are parametrized by the mean value and the standard deviation. Log-normal distributions are considered to avoid physically impossible negative realizations. We employ standard random number generators [36] to obtain the sequence of samples required for Monte Carlo sampling.
In the case of heterogeneous random fields,m i (x), we employ stationary log-normal random fields whose spatial correlation is represented by the auto-correlation function where x 1 and x 2 are two points in a background domain which is larger than all fracture length realizations, and l m i is the correlation length. To generate samples of the random field,m i (x), it must be discretized. To obtain a discretization the log-normal random field is considered as the exponential of an underlying stationary normal random fieldg i : The statistical moments of the underlying Gaussian distribution can be expressed in terms of those of the random fieldm i (x) by and the auto-correlation function as Discretization of the underlying Gaussian random fieldg i (x) is then achieved by means of the Karhunen-Loève expansiong wherez is a vector of n independent standard normal random variables, and where ξ j and r j (x) are the eigenvalues and eigenfunctions corresponding to the spatial covariance function σ 2 g i ρ g i (x 1 , x 2 ), respectively. We discretize the eigenfunctions in space by means of a uniform linear finite element discretization over the background domain, which results in a generalized eigenvalue problem that we solve using a direct eigenvalue solver (see Appendix 2 for details). In Figure 4a the first 12 numerically determined eigenfunctions are shown for the auto-correlation function (25) with l m i = 10.
The log-normal random fieldm i (x) is then obtained by back-substitution of (29) into the transformation (26) to yieldm Realizations of the random fieldm i (x) can now be generated by sampling a sequence of n independent standard normal random variables. In Figure 4b we show the approximate auto-correlation function for various numbers of Karhunen-Loève terms, which conveys that the approximate autocorrelation function converges to (28) by increasing the number of random variables n. Note that the number of random variables required to obtain a fixed accuracy increases as the correlation length decreases.

Numerical simulations
In this section we present numerical results based on the methodology presented above. In Section 5.1 we validate our methodology by consideration of the benchmark results presented in the comparative study by Warpinski et al. [37]. In this section we demonstrate the necessity to use a tip enrichment function and enforcement of volume conservation, and we study the influence of the mesh size and time-step size on the numerical results. In Section 5.2 the sensitivity of the observables -in particular the fracture length and aperture -to the uncertain model parameters are studied, which serves as a starting point for the stochastic simulations discussed in Section 5.3.
In the stochastic setting the uncertain model parameters are represented by discrete random fields.

Deterministic benchmark
To demonstrate the validity of the presented methodology we consider the benchmark case studied by Warpinski et al. [37], which is based on a staged-field experiment of the Gas Research Institute [37, p. 26]. The considered model parameters are assembled in Table 1. The injection rate is gradually increased until the indicated value, and then held constant for 200 minutes. The material parameters resemble that of a tight gas sand reservoir, for which spurt losses are omitted. In Figure 5a we show the evolution of the fracture in time, where it should be noted that the height of the fracture, H, is constant. Since the width is symmetric with x axis, in the figure , just the half width aperture is plotted. The shown results are based on a mesh size of ∆x = 1 m and a time step size of ∆t = 1 s. As we will study in detail below, these results are objective with respect to these numerical parameters. Figure 5b shows the corresponding increase in fracture length and fracture mouth opening over time. The observed time evolution corresponds well with the results reported for various simulators in [37]. It is noted that the reported results in [37] vary significantly as a result of variations in model assumptions and simulation frameworks. The fracture length of 1429 m as computed here also corresponds reasonably well with the analytical model in [25], which predicts a fracture length of 1340 m. Note that in the absence of leakoff our model predict a fracture length of 1730 m. This stipulates that leak-off is appropriately represented in our numerical simulations. The fracture length and fracture opening computed by our methodology are in the higher part of the spectrum of simulators considered in [37] and analytical models, which we attribute to the enforcement of the volume constraint, which will be discussed in detail below.
The above benchmark results are based on our methodology including the enrichment of the tip functions and the enforcement of the global volume conservation constraint. We have found both aspects to be essential in order to establish numerical results with an acceptable level of accuracy for meshes and time step sizes that are computationally tractable in the scope of stochastic simulations. In Figure 6 we study the behavior of the global volume conservation without and with Lagrange multiplier constraint, respectively. The total volume rate -which is the sum of the leak-off rate and fracture-widening rate -should equate to the input flow rate. Note that in the absence of the Lagrange multiplier constraint, a significant mismatch between the total rate and the inflow rate is observed. The presented figure is based on on a mesh size of ∆x = 1 m and a time step size of ∆t = 1 s . The mismatch between the rates depends on these discretization parameters, as it originates from significant errors in the local volume balance in the finite element discretization (15). These local inaccuracies in the finite element solution are closely related to the highly nonlinear character of the constitutive relation. By enforcing global conservation of volume using a Lagrange multiplier -as shown in Figure 6b -the global loss of volume is rigorously resolved. As observed the total volume rate in this case matches that of the inflow rate up to a specified numerical tolerance.  (12) is used to compute the tip speed. In the absence of enrichment the tip speed cannot be obtained by this equation, as the adequate singular tip behavior is then not represented in the discrete solution space. The speed results presented in Figure 7a are based on the finite difference approximatioṅ Note that both equation (12) and this finite difference approximation are evaluated after the subiteration procedure has converged. From 7b it is observed that enriching the discrete solution space benefits the simulation, as the computed tip velocity is very smooth in comparison to that computed using a solution space without enrichment. The enforcement of the volume constraint by means of a Lagrange multiplier, and the representation of the tip behavior by means of an enrichment function, result in solutions with a level of accuracy that enable studying stochastic variations. In Figures 8,9 we show the dependence of the results on independent variations in the time step size and mesh size, respectively. Note that in all simulations we consider a duration of 50 s only, in order to make the converge studies feasible in terms of computational effort.
From Figure 8 the simulation for a mesh size of ∆x = 1 m are studied using three time step size, viz. ∆t = 1.0, 0.5, 0.25 s. From both the length evolution plot and the tip speed evolution plot it is observed that the variations with the time step size are very limited. The most notable difference is observed at the onset of fracture, where the maximum tip speed for ∆t = 0.25 s is observed to be 2% higher than that using ∆t = 1.0 s. This difference is significantly smaller once steady propagation occurs, e.g. at t = 50 s, where the difference is only 0.8%. Since the length of the fracture is generally not dominated by the onset phase, the observed variation in fracture length is generally also very small. At t = 50 s, the fracture length for ∆t = 0.25 overestimates that of ∆t = 1.0 by 0.5%. Although not presented here for the sake of brevity, similar results can be established for other indicators such as the fracture mouth opening.
The dependence of the fracture length and tip speed is depicted in Figure 9. From the tip speed evolution it is observed that although a mesh size of ∆x = 2.0 m correctly mimics the tip speed behavior, fluctuations in the speed can be observed as a consequence of the mesh coarseness. These fluctuations can be attributed to the fact that due to the limited number of elements in this simulation (i.e., only 10 elements at t = 50 s) the spatial discretization errors resulting from the employed moving mesh approach are significant. Upon mesh refinement these effects are observed to vanish. The maximum tip speed at the onset of fracture for a mesh with ∆x = 0.5 m is observed to be only 0.5% lower than that using a twice as coarse mesh. When the fracture is steady at t = 50 s, this relative difference is observed to be even smaller. In terms of the fracture length at t = 50 s, the mesh with ∆x = 0.5 m overestimates that of the mesh with ∆x = 1.0 m by a mere 0.5%. In the context of the sensitivity studies and stochastic simulations considered in the remainder of this section, it is essential that the numerical errors do not pollute the results. This means that the variations discussed above should be negligible in comparison to the stochastic variations in the input parameters. On the other hand, making the mesh sizes and time step sizes too small will dramatically increase the computational effort. Herein we opt to use a mesh size of ∆x = 1.0 m and a time step size ∆t = 1.0 s, with which we achieve a good balance between numerical accuracy (see above) and computational effort.

Sensitivity analysis
In this section we leverage the deterministic model outlined above to identify the input factors that drive the variation in the output. As output observables we consider the fracture length and fracture aperture. As input parameters we consider the model parameters that cannot be established with a high degree of certainty in reality, viz. the plane strain modulus E , the fracture height, H, and the leak-off coefficient, c l . In our sensitivity analysis we independently vary the input factors and study their impact on the output observables. This screening exercise is often considered as the first step in a forward uncertainty analysis, since it identifies the dominant sources of randomness.
In Figure 10 we consider the effect of the plane strain modulus on the fracture geometry, while keeping all other model parameters unchanged. A range of plane strain moduli from 1 × 10 3 MPa to 1 × 10 4 MPa is considered. From Figure 10 is observed that a stiffer formation results in a longer and narrower fracture, compared to a more compliant formation. Considering the nature of the model -which revolves around the conservation of volume -this makes sense, since in the case of a stiff formation fracture propagation is favored over fracture widening. From Figure 10a we observe that the dependence of the output observables on the plane strain modulus is highly nonlinear, in the sense that the rate of change of the output observables is significantly lower than that of the plane strain modulus. For example, an increase of 45% in the fracture length is observed when increasing the plane strain modulus by a factor of 8. This same increase in plane strain modulus is observed to reduce the fracture mouth opening by only 35%. As a result, the well pressure -which is proportional to the product of the fracture mouth opening and plane strain modulus -increases with an increase in formation stiffness, which is inline with the experimental results [37]. Another observation in this sensitivity analysis is that the response of the observables is non-symmetric, in the sense that the rate of chance of the length and fracture mouth opening for stiffer formations is smaller than for compliant formations. Figures 11 shows the dependence of the output observables on a range of leak-off coefficients, ranging from the impermeable case (c l = 0 m/s 1/2 ) to c l = 5e −5 m/s 1/2 , which is 50% more than the value taken in GRI experiment [37] which is based on tight gas sands. As anticipated from the conservation of volume, an increasing leak-off coefficient yields a shorter and narrower crack. Increasing the deterministic value considered in the previous section by a factor of 5 yields a decrease in fracture length of 14% and a decrease in fracture mouth opening of 7%. From Figure 11b it is observed that the fracture profile shape is virtually insensitive to the leak-off coefficient. In contrast to the dependence on the plane strain modulus considered above, the rate of change of the observables is practically constant for the considered range of leak-off coefficients. In Figures 12 the variation of the output observables for fracture heights ranging from 25 m to 95 m is observed. Doubling the fracture height from 25 m reduces the fracture mouth opening by 15% and the fracture length by 44%. As expected from volume conservation (where in this particular case leak-off effects are not pronounced) the product of these two observables is approximately reduced by a factor of two. The direct impact of the fracture height on the volume conservation model results in a strong sensitivity of the output observables. The non-symmetry of the response observables is consistent with the expected behavior in the extreme cases, for which a extremely large height should yield a very short and narrow crack, and a extremely small height (which is a violation of the model assumptions) should yield an extremely long and wide fracture.

Stochastic setting
In this section we present the results of the Monte Carlo simulations. In Section 5.3.1 we first study the stochastic results for the case where each of the uncertain input parameters is varied independently, which closely connects this section to the sensitivity analysis presented above. The stochastic analysis presented here provides additional insights into the evolution of the randomness  Table 2: Basic statistical properties for variation of E for different times in time, and on the dependence of the uncertain observables on the magnitude of the random input variables. In Section 5.3.2 we then consider a case of a random field for the plane strain modulus, which provides additional insights into the dependence of the observables on the spatial correlation of the uncertain input parameter. The reported sample sizes are all based on the estimate (24) with a confidence level of 95% for the estimator of the mean fracture length.

Independent variation of uncertain parameters
We first consider the plane strain modulus, E , to be the only uncertain parameter. Since the plane strain modulus is positive by definition, we represent this uncertain parameter by a log-normal distribution with mean value µ E = 6.13 × 10 4 MPa and coefficient of variation V E = 15% . In Figure 13 we show the evolution of the fracture length in time using a sample size of N = 304, where the mean value is represented by the solid line, and the shaded area indicate the 98% confidence interval of the fracture length (V L ≈ 0.1). Typical distributions for the fracture length and fracture mouth opening at t = 100 s are displayed in Figure 14, from which it is observed that the distributions are unimodal. The error bars in Figure 13 show that the standard deviation of the fracture length (and with that its confidence interval) increase proportionally with the mean. This observation is confirmed in Table 2, from which it is observed that the coefficient of variation of the length and fracture mouth opening only varies moderately.
In Figure 15 we study the dependence of the results on the coefficient of variation of the plane strain modulus at two time instances. Note that the used sample size is different for each coefficient of variation in order to achieve the same confidence level for the mean estimator of the fracture length. At both time instances we observe the coefficient of variation of the fracture length to be proportional to that of the input parameter. This behavior can be explained by considering the dependence of the mean and the standard deviation on the coefficient of variation of the input, as displayed in Figure 16. From this figure it is observed that the standard deviation of the fracture length increases with an increasing variation of the plane strain modulus. Using the results from the sensitivity analysis presented above we verify that the observed behavior is in good correspondence with that expected from first-order perturbation analysis [38]: The observed decrease in the mean value is also observed to be in good correspondence with the results of the perturbation analysis, as shown in the Figure 15a, where the second-order expectation of mean is computed using [38] µ L ≈ L| E =µ E + 1 2 For the case considered here the variation in the mean of the fracture length is moderate, as a result of which the observed behavior of the coefficient of variation in Figure 15a is governed by the behavior of the standard deviation of the fracture length. Similarly, the results for the fracture mouth opening are also observed to correspond well with the sensitivity analysis.
In Figure 17a we show the evolution in time of the fracture length for the case in which the leak-off coefficient is described by a log-normal distribution with mean value µ c l = 9.84 × 10 −6 m/s 1/2 and coefficient of variation V c l = 50%. The probability distributions at t = 100, s are shown in Figure 18. In Figure 19 we study the relation between the coefficients of variation of the observables and that of the leak-off coefficient at two time instances, where the sample sizes have been selected in accordance with the confidence level of the mean estimator of the fracture length. It is observed that the coefficients of variation of the output observables are not changing significantly in time. The observed relation between the coefficients of variation of the input and output is explained by the fact that the mean value is effected minimally by the coefficient of variation of the leak-off coefficient, while the standard deviation increases proportionally with it. From Figure 19 it is observed that the behavior of the mean and standard deviations of the observables is in good agreement (with a maximum of 3% variation which can be attributed to the approximations made while computing perturbation results, see Appendix 4)with the perturbation results following from the sensitivity analysis. We finally consider the independent random variation of the fracture height, which is represented by a log-normal distribution with mean µ H = 51.8 m and coefficient of variation V H = 50%. A sample size of N = 553 (with a confidence level of 95% ) is selected to compute the time evolution of the fracture length as shown in Figure 17b, for which two of the probability distributions at time t = 100 s are shown in Figure 21. The coefficient of variation is observed not to be subject to significant changes in time, which is confirmed by comparison of the relation between the input  Figure 22. From Figure 23 we observed that both the mean and the standard deviation of the observables match well with the results from the sensitivity analysis.
Comparing the effects of the random input variable on the fracture length reveals that it is more sensitive to randomness in the plane strain modulus than to randomness in the leak-off coefficient, in the sense that a coefficient of variation in the fracture length of V L ≈ 2.5% corresponds to a coefficient of variation of the plane strain modulus of V E ≈ 15% and a coefficient of variation of the leak-off coefficient of V c l = 50%. The sensitivity to the fracture height is observed to be even strong, in the sense that V L ≈ 10% corresponds to coeffients of variation of V E ≈ 50% and V h f = 15% for the plane strain modulus and fracture height, respectively.

Heterogeneous random plane strain modulus field
We now consider the case in which the plane strain modulus is described by a random field, E (x), with mean µ E = 6.13 × 10 3 MPa and coefficient of variation V E = 50%. We vary the length scale of the auto-correlation function Equation (28) from l E = 5 m to l E = 25 m for t = 100 s. The random field for the plane strain modulus is discretized using 12 Karhunen-Loève modes, which is sufficient for the representation of the random field corresponding to the smallest correlation length considered. In Table 3 the statistical moments of the observables are collected based on a Monte-Carlo simulation with N = 384, which is in accordance with a 95% confidence level for the mean estimator in the fracture length.
The most notable observation from the results in Table 3 is that the coefficient of variation of the output observables is significantly higher than in the case of a homogeneous plane strain modulus with equal coefficient of variation (see Table 3). To better understand this observation, in Figure 24 we perform a closer inspection of the realizations that lead to Table 3. In the rows of this figure we collect the Monte-Carlo results for the correlation lengths reported in Table 3, starting with the smallest correlation length on top. In the second and third column we show the probability distributions for the fracture length and the fracture mouth opening. In the first column we show the plane strain modulus field that leads to three distinct realizations in the sample, viz. the smallest fracture length, the largest fracture length, and the fracture length closest to the mean value.
We observe that the realizations of the plane strain modulus field that lead to the smallest fracture lengths in all cases correspond to the situation in which the elastic modulus is very small near the well. When this happens the injected fluid causes fracture widening near the well, rather than fracture propagation into the formation. More generally, in the case of heterogeneous fields, local zones in which the formation is very compliant can lead to blockage of propagation, as the injected fluid volume can be locally accumulated in this zone. Long fractures are obtained in the case that the plane strain modulus is large near the well, and high (in a spatially averaged sense) compared to the mean value. In such situations the blockage of propagation due to a compliant zone does not occur.
In terms of the dependence of the results on the correlation length it is observed that the mean fracture length decreases as the correlation length decreases. This is explained by the fact that in the case of a smaller correlation length, the chance of a locally compliant zone in the formation increases. The blockage of propagation in such zones is then more frequent, which leads to a reduction in fracture length expectation. From Table 3 we moreover observe a moderate increase in coefficient of variation of the fracture length as the correlation length increases.
From the distributions of the fracture mouth opening in Figures 24c1-24c5 we observe a notable difference in comparison to those for the homogeneous random plane strain modulus case (Figure 14b). In the homogeneous case there exists a strong correlation between the fracture length and the fracture mouth opening, in the sense that long cracks are narrow by virtue of the fact that their volume is similar (assuming leak-off effects to be limited). Figure 10 in the sensitivity study clearly confirms this observation. Although the fracture length and fracture width in the case    of a heterogeneous field are not uncorrelated, the fracture mouth opening is strongly influenced by the local plane strain modulus near the well. Since the fracture opening in the PKN model depends locally on the plane strain modulus, the log-normal distribution of the plane strain modulus reflects directly on that of the fracture mouth opening, as can be seen in the third column of Figures 24c1-24c5. The sensitivity of the fracture mouth opening to local variations in the plane strain modulus field also results in coefficients of variation that are significantly higher than those in the homogeneous case.

Conclusions
We have presented a sampling-based stochastic analysis of the hydraulic fracturing process based on the Perkins-Kern-Nordgren (PKN) model. The consideration of this model is motivated by the fact that in the deterministic case high-accuracy solutions can be computed with limited computational effort, which makes its application in the context of direct Monte Carlo sampling practical. Although this model significantly simplifies the hydraulic fracturing process, it bears practical relevance, especially for fractures in the viscosity-dominated regime. A limitation of the model of particular relevance in this study results from the local elasticity relation in the PKN model, which restricts its application to low-frequency spatial variations of the model parameters.
In order to compute high-fidelity solutions that do not pollute stochastic analyses with numerical errors, a moving-mesh finite element method is developed. The employed backward-Euler time integration scheme is supplemented with a sub-iteration technique, such that the mesh propagation relation becomes implicit. The non-linearity of the model is solved using Newton iterations.
We have performed detailed mesh size and time integration step convergence studies. We have found that in order to attain solutions with acceptable accuracy in the context of the stochastic  Figure 26: (a) )µ x f and µ wmax with different variation of E at t = 100 s for l E = 20 m, and (b) σ x f and σ wmax with different variation E at t = 100 s for l E = 20 m analysis considered herein, the finite element methodology had to be enhanced. First, the global conservation of volume was found to be significantly violated due to the highly non-linear character of the model. The observed loss of volume led to significant underestimation of the fracture length. To circumvent this problem, volume conservation was enforced by means of a Lagrange multiplier approach. Second, the weakly singular behavior of the fracture opening and pressure at the tip was found to be troublesome in the case of a standard finite element basis. On one hand the improper representation of the singularity by the basis required the use of an ad hoc tip velocity relation. On the other hand, the mesh resolution of the uniform finite element mesh was found to be insufficient. These issues were resolved by enrichment of the standard finite element space with a singular tip function.
Our finite element simulations show very good agreement with results reported in literature for a realistic test case. The sensitivity of the fracture evolution process with respect to various random input parameters was studied. From the conducted direct Monte Carlo simulations it was found that the mean and standard deviation of the fracture length and fracture mouth opening correspond well to those values obtained using perturbation theory. This observation conveys that -at least for the test case considered -linearization of the model provides meaningful information on the behavior of the stochastic moments, despite the complexity of the model and its solution procedure.
To demonstrate the suitability of the developed methodology for studying random heterogeneities in formation properties we have considered a test case in which the formation stiffness was described by a random field. The random dimension was discretized using a Karhunen-Loéve expansion, while the spatial dimension was discretized using a linear finite element mesh on a background domain. The sampling results demonstrate that the response uncertainty is amplified by the heterogeneous character of the random material property field. For the fracture length this is explained by the fact that fracture propagation is sensitive to local variations in the elastic properties of the formation. For the fracture mouth opening an even stronger amplification is observed as a consequence of the fact that the fracture opening is directly related to the material property. Although this observation can be explained well based on the structure of the PKN model, it requires further study to understand to what extend a similar conclusion can be drawn for more rich hydraulic fracturing models.
Although the results presented herein provide excellent insight into the primary characteristics of the stochastic behavior of the hydraulic fracturing process, it is evident that more detailed information can be obtained by a more versatile model and simulation strategy. In particular the PKN model does not rely on a fracture mechanics model based on the material's fracture toughness, which restricts the scope of this work to fractures in the viscosity-dominated regime.
When considering uncertainty quantification using physically richer models it will remain key to not pollute the results with numerical errors, which will inevitably lead to computationally very expensive Monte Carlo methods. The use of alternative stochastic techniques, such as e.g. the perturbation method can be expected to yield meaningful results at a much lower cost than direct Monte Carlo sampling.

A Karhonen-Loeve Expansion and its discretizion using finite elements
Here, we discuss the discretization of Karhunen-Loève expansion in a finite element setting. The Karhunen Loéve expansion of a normal random field is expressed as: wherez i are stochastic dependence or uncorrelated standard normal variables and g i (x) are the eigen functions of the following integral equation [38].
where EE (x 1 , x 2 ) are covariance kernal of the deterministic function of the random fieldG(x). Furthermore, the random process can be expressed as a direct sum of orthogonal projections of the basis functions g i (x) which are proportional to the corresponding eigenvalues.
Thus, it has a spectral decomposition The eigenvectors and eigenfunctions can be solved using finite elements. The eigenfunctions are represented in terms of the shape functions. The infinite number of eigenfunctions are approximated to finite number and are obtained by expanding in finite elements and transforming the eigenvalue problem into an algebraic eigenvalue problem [39].
Discretization the integral equation we get, Galerkin formulation results in generalized matrix eigen value problem of the form Thus, random field in finite elements for a normal field can be written as B Length and maximum aperture of the fracture at different times We consider the benchmark case studied by Warpinski et al. [37], which is based on a stagedfield experiment of the Gas Research Institute [37, p. 26]. The considered model parameters are assembled in Table 1. The injection rate is gradually increased until the indicated value, and then held constant for 200 minutes. The material parameters resemble that of a tight gas sand reservoir, for which spurt losses are omitted. The shown results are based on a mesh size of ∆x = 1 m and a time step size of ∆t = 1 s and rounded off to 4 significant digits.

C Perturbation results
Using the results from the sensitivity analysis from presented above we verify that the observed behavior is in good correspondence with that expected from first-order perturbation analysis [40]: The observed decrease in the mean value is also observed to be in good correspondence with the results of the perturbation analysis, as shown in the Figure 16a, where the second-order expectation of mean is computed using the following equation (43) µ L ≈ L| E =µE + 1 2 For all the calulations, we compute the first dervatives and second dervatives using finite differences scheme.
In case of E , the corresponding values of length and width were from 5 × 10 10 P a and 7 × 10 10 P a. We get the following In case of cl, the corresponding values of length and width were from 9 × 10 −6 m/s 1/2 and 1 × 10 −5 m/s 1/2 . We get the following In case of H, the corresponding values of length and width were from 45 m and 60 m. We get the following