Handling spatial heterogeneity in reservoir parameters using POD

Subscriber access provided by UNIV OF DURHAM

Article

Handling spatial heterogeneity in reservoir parameters using PODbased EnKF for model-based feedback control of hydraulic fracturing Abhinav Narasingam, Prashanth Siddhamshetty, and Joseph Sang-Il Kwon Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b04927 • Publication Date (Web): 26 Feb 2018 Downloaded from http://pubs.acs.org on March 6, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Handling spatial heterogeneity in reservoir parameters using POD-based EnKF for model-based feedback control of hydraulic fracturing Abhinav Narasingam, Prashanth Siddhamshetty, and Joseph Sang-Il Kwon∗ Artie McFerrin Department of Chemical Engineering, Texas A&M University, College Station, TX Texas A&M Energy Institute, Texas A&M University, College Station, TX E-mail: [email protected]

Abstract Accurate characterization of reservoir properties is of central importance to achieve a desired fracture geometry during a hydraulic fracturing process. However, the estimation of spatially varying geological properties in hydraulic fracturing is inherently ill-posed due to a limited number of measurements. In this work, parameterization is performed to reduce the dimensionality of spatially varying Young’s modulus profiles via proper orthogonal decomposition (POD), and a data assimilation technique called Ensemble Kalman filter (EnKF) is used to estimate the parameter values in the reduced low-dimensional subspace. Through a series of simulation results, it is demonstrated that the POD-based EnKF technique provides a process model with updated spatially varying geological parameters, which is able to make an accurate prediction of the fracture propagation dynamics in hydraulic fracturing. Next, we use the updated high-fidelity process model in a model predictive control framework to construct a closed-loop system under the proposed model-based feedback control system that achieves uniform proppant concentration at the end of pumping in a hydraulic fracturing process. ∗ To

whom correspondence should be addressed

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction In recent years, the extraction of underground resources such as shale gas and oil, which are trapped in ultra-low permeability formations, has become economically viable due to the advances in wellstimulation techniques such as hydraulic fracturing. 1 In hydraulic fracturing, the first stage is called perforation, in which a well is drilled and a wire equipped with explosive charges is dropped into the well to create initial fracture paths. Next, a high pressure fluid mixture consisting of water and viscosifying agents such as proppant (e.g., sand or ceramic materials) is introduced into the wellbore at a sufficiently high flow rate to break the rock and propagate the fractures in the rock formation. In the final stage, the pumping is halted and the remaining fluid seeps into the reservoir, trapping the proppant between the fracture walls. This process will create a highly-conductive channel that facilitates the effective transport of oil and gas to the wellbore. 2 The final fracture conductivity depends on the propped fracture geometry and the proppant distribution throughout the fracture. Recently, several efforts have been made for the numerical analysis of hydraulic fracturing and the design of optimal control strategies to achieve the desired fracture geometry and uniform proppant concentration at the end of pumping. 3–7 The above-mentioned work employed reservoir models where the rock mechanical properties were assumed to be known and spatially homogeneous. However, field data indicates that even within the same rock formation, the performance of hydraulic fracturing can be notably different. 8,9 This variability can, in part, be attributed to the spatially varying rock mechanical properties such as the Young’s modulus. As a consequence, one of the key tasks, which has to be performed prior to the model-based controller design, is the characterization of spatially varying rock mechanical properties. This requires the availability of process data. Currently, in the field, a small-scale preliminary experiment called the mini-frac test is performed to collect preliminary data required for pre and post-closure fracture analysis. This data can further be used for parameter estimation to characterize the variability in the spatially varying geological properties. A well-known technique for estimation of states and parameters of a nonlinear system is the 2


Page 2 of 43

Page 3 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ensemble Kalman filter (EnKF). The EnKF was introduced initially by Evensen, 10 and has been used successfully in several applications. 11–14 In summary, the EnKF can be interpreted as a sequential Monte-Carlo approach initiated by choosing an ensemble of model parameters sampled from an initial distribution generated using a priori (geostatistical) assumptions. Then, the ensemble propagates through a high-fidelity model and is updated using the available measurements. After all the measurements are assimilated, the mean of the final ensemble gives the approximate representation of the unknown spatially varying parameter profile. In a typical hydraulic fracturing process, real time measurements may include downhole pressure and micro-seismic data which will be processed to provide the fracture width at the wellbore and the fracture length, respectively. 15–17 Microseismic monitoring (MSM) is the detection of signals produced by small sesimic events and is used for determining the location and geometry of fractures. When a crack occurs in a rock formation, energy is released, referred to as microseismicity, and some of the energy travels away from the source through the surrounding rock as seismic waves. These waves temporarily deform the material as they travel in the form of P-waves, or primary waves, which are relatively fastpropagating waves and S-waves, or secondary waves, which are slower-propagating shear waves. MSM operations incorporate arrays of three-component geophones or accelerometers deployed in a nearby monitoring well that pick up these seismic events. The distance of the originating source of a microseismic event to each geophone is calculated based on the difference in arrival time between P- and S- waves and a previously calibrated seismic-wave velocity model. The detection results of the same event by three or more geophones allow determination of the location of the event source in three-dimensional (3D) space. Analysts interpret these MSM locations to show induced fracture length, height and azimuth. Downhole pressure data can be used to estimate the fracture width near the wellbore via pressure-width relationship equations. The hydraulic fracturing process is described by a system of nonlinear partial differential equations with an added complexity that the spatial domain varies with time as the fracture propagates. Therefore, a large number of discrete elements are required to provide a detailed description of the reservoir. This results in the formulation of a large-scale inverse problem where the number

3



of parameters to be estimated far exceeds the number of available field measurements, resulting in an ill-posed problem. To deal with this unidentifiability issue, the original problem must be reformulated and one way to achieve this is to reduce the number of unknown model parameters by parameterization. Additionally, it is of critical importance to preserve the spatial features in the geological properties during the parameterization process. Proper orthogonal decomposition (POD) 18,19 is a common transformation-based method that provides a way to reduce the number of parameters while still maintaining the key spatial features. Although it has been used extensively for reservoir characterization in history-matching problems, 20–23 its application to estimation of spatially varying geological properties in hydraulic fracturing has not received much attention. In this work, we use an integrated framework to parameterize the unknown spatially varying Young’s modulus profile via POD, keeping the key spatial pattern in the geological properties intact, and to update the statistical information using the available process measurements via EnKF. The high-fidelity model of hydraulic fracturing with the updated spatially varying Young’s modulus profile is simulated to generate the data required for model order-reduction. The local dynamic mode decomposition with control (LDMDc) methodology, which is developed in our previous work, 24 is then applied to the high-fidelity simulation data to derive reduced-order models. It has been demonstrated that the proposed LDMDc methodology provides accurate approximations to the high-fidelity solutions of nonlinear distributed parameter systems, particularly in the presence of spatial domains that evolve with time. These reduced-order models are used in the model predictive controller (MPC) to drive the fracture geometry and the proppant concentration at the end of pumping to their respective set-point values, which are determined by the rock permeability and the amount of proppant to be injected, by manipulating the flow rate and inlet proppant concentration at the wellbore. The rest of the paper is organized as follows: initially, the high-fidelity model of hydraulic fracturing processes is introduced to show how spatially varying Young’s modulus profile is taken into account in the model. Then, a description of the POD-based EnKF methodology for the parameterization and estimation of the spatially varying Young’s modulus profile is presented. This is

4


Page 4 of 43

Page 5 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


followed by the section on the development of a model-based feedback control system to compute the optimal pumping schedule. Finally, the closed-loop simulation results of hydraulic fracturing are presented to demonstrate the significance of incorporating spatially varying parameter profiles in the feedback control system.

Dynamic modeling of Hydraulic Fracturing Hydraulic fracturing can be modeled as a dynamic process describing the fracture propagation, fluid flow, and proppant transport inside the fracture. The fracture propagation is modeled based on the following assumptions: (1) fracture propagation is described using the Perkins, Kern, and Nordgren (PKN) model 25 ; (2) the formation layers above and below are where the fractures have sufficiently large stresses such that the vertical fracture is confined within a single horizontal rock layer, i.e, the height of the fracture is constant in the vertical direction; and (3) because the fracture length is much greater than its width, the fluid pressure across the vertical direction is constant. A brief description of the model equations is presented below. For a detailed description we refer the readers to Ref. 26,27 By simultaneously taking into account the conservation of local fluid mass, lubrication theory, elasticity equation, and a model describing the rate of fracture fluid leak-off into the formation, a nonlinear partial differential equation (PDE) with a time-dependent spatial domain is formulated as follows: 2 2Cleak ∂ A ∂ qz + +Hp =0 ∂t ∂z t − τ(z)

(1)

where A = πW H/4 is the cross-sectional area of the elliptic fracture, H is the fracture height, z ∈ [0, L(t)] is the time-dependent spatial coordinate in the lateral direction, W (z,t) is the fracture width, Cleak is the overall leak-off coefficient, t is the elapsed time since the hydraulic fracturing process was initiated, τ(z) is the time at which the fracture has arrived at the position z for the first time, and qz (z,t) is the local flow rate in the horizontal direction, which is related to the fracture

5



Page 6 of 43

width using lubrication theory and elasticity equation as follows: π 3 ∂W 4 dE(z) qz = − + E(z)W W 128η(1 − ν 2 ) dz ∂z

(2)

where E(z) is the Young’s modulus of the formation which varies with the spatial coordinate z, η is the fracture fluid viscosity, and ν is the Poisson ratio of the formation. The spatially varying Young’s modulus, E(z), describes the in-situ stress environment within the rock formation and has a direct impact on the geometry of the induced fracture. Therefore, it is important to accurately characterize it. Remark 1 Cleak is a measure of the flow resistance of the fluid leaking off into the formation during hydraulic fracturing. 28 It defines the three types of linear flow mechanisms encountered with fracturing fluids. These are: (1) viscosity and relative-permeability effect; (2) reservoir-fluid viscosity-compressibility effects; and (3) wallbuilding effects. The first two mechanisms involve coefficients which can be calculated from reservoir data and fracturing-fluid viscosity with the aid of presently available formulae. The third case involves fluid-loss coefficients for fluid-loss additives, which must be experimentally determined. In summary, the leak-off coefficient considered in this work is a constant whose value depends on the specific type of reservoirs and fracturing-fluids used in the operation. At the wellbore z = 0, the injected flow rate qz is specified, and at the fracture tip z = L(t), the fracture is always closed (i.e., the fracture width is zero), leading to the following two boundary conditions: qz (0,t) = Q0

and W (L(t),t) = 0

(3)

where Q0 is the fluid injection rate at the wellbore (i.e., the manipulated input). Initially, the fracture is closed, leading to the following initial condition:

W (z, 0) = 0

6


(4)

Page 7 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


By considering simultaneous advection in the lateral direction and settling in the vertical direction, proppant transport is modeled as follows: ∂ (WC) ∂ + (WCVp ) = 0 ∂t ∂z

(5)

C(0,t) = C0 (t) and C(z, 0) = 0 where C(z,t) is the proppant concentration inside the fracture and C0 (t) is the injected proppant concentration at the wellbore. The velocity of an individual proppant particle, Vp , is given by the summation of the fluid velocity, V , and the gravitational settling velocity, Vs , as follows: 29,30

Vp = V − (1 −C)Vs

Vs =

(1 −C)2 (ρsd − ρ f )gD2 101.82C 18η

(6)

(7)

where ρsd is the proppant particle density, ρ f is the pure fluid density, g is the gravitational acceleration constant, and D is the proppant diameter. The fracture fluid viscosity is related to the proppant concentraion via the following empirical expression: 31

η(C) = η0

C 1− Cmax

!−α (8)

where η0 is the pure fluid viscosity, α is an exponent in the range of 1.2 to 1.8, and Cmax is the maximum theoretical concentration determined by Cmax = (1−φ )ρsd where φ is the proppant bank porosity. The evolution of proppant bank height, δ , by the settling flux is described by, 32 d(δW ) CVsW = dt (1 − φ )

(9)

where the initial condition is that δ (z, 0) = 0. Please note that in the above equations we have not considered the effect of fracture wall roughness on the fluid transport. There have been several studies examining these effects of frac7



ture aperture and fracture roughness on the mechanical (elastic) and hydraulic (fluid transport) response of a fractured rock. For example, the effect of roughness on mechanical response includes the stress dependence, which has a significant impact on seismic monitoring, 33,34 that is shown to be useful in potentially estimating the fracture size/length. 35,36 Additionally, as the contact area between the fracture surfaces increases with the normal stress, or as the fracture surfaces become rough, the fluid flow pattern also changes affecting the hydraulic response of the fractured rock. Under the rough wall conditions, the fluid will not flow through the entire fracture surface but rather through a few preferred paths (channels). Several researchers have come up with models 37–42 for studying the characteristics of fluid flow within rough fracture walls. All these models propose modifications to the lubrication (Boussinesq) equation which relates the fracture aperture and the pressure gradient via a cubic law. Ideally, under the assumption of smooth fracture walls, the fracture aperture in the equation is considered constant but in the presence of roughness, the proposed models treat the fracture aperture as a function of fracture roughness using empirical correlations. However, since all the modifications proposed will take effect only in modeling the relation between the fracture aperture and the fracture roughness, it does not significantly affect the data assimilation process adopted later in this manuscript. The procedure would remain the same with some changes to the above dynamic model of the hydraulic fracturing process. The system of Eqs. (1)-(9) is solved using a novel in-house high-order discretization scheme previously developed in Ref. 27 Fig. 1 shows the evolution of the fracture width and proppant concentration for Q0 = 0.03m3 /s starting from the initial condition. It can be observed that the growth rate of the fracture width is very high in the beginning but it slows down with time. At the same time, the fracture length grows steadily at a constant rate. It is important to note that the zero values in the above figure indicate the fracture has not yet propagated to that location. After discretization in space and time, each element contains two states, namely the fracture width and the proppant concentration. Since the process is characterized by a set of highly coupled nonlinear equations defined over the time-dependent spatial domain, a large number of discretized elements are necessary to describe the hydraulic fracturing process accurately. This leads to a large-scale

8


Page 8 of 43

Page 9 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


system in both the state and parameter space. As a result, reduction in the number of parameters via parametrization is a crucial step in reservoir characterization for developing accurate reducedorder dynamic models 43 that can be used for designing state estimators and model-based feedback control systems. 24

POD-based EnKF This section provides a brief description of particular paramterization and data assimilation techniques essential for characterizing the spatially varying Young’s modulus profile of the reservoir. In this work, we employ the EnKF, which is a statistical Monte-Carlo simulation where the ensemble of process model states are evlolved forward in time with the ensemble mean as the best state estimate 44 and the ensemble spread as the error variance. It combines a stochastic model with a prior assumption about the states, process noise, and measurement noise to derive the first two statistical moments of the states’ posterior distribution. In other words, it is a sequential state estimation method for dynamical systems. It has also been widely used to estimate process model parameters. 45–47 In this work, we use it for solving the inverse problem to estimate the spatially varying Young’s modulus profile of the reservoir. In a practical implementation, the EnKF involves a forecast step which predicts the mean and covariance of the process model parameters at the next time step using the high-fidelity process model of hydraulic fracturing described by Eqs. (1)-(9). While the true values of the spatially varying parameters are unknown, we can assume that we have access to the first two statistical moments (i.e., mean and covariance) of the underlying geological features of the formation based on historical data. Please note, in reservoir characterization problems, the input reservoir properties are unknown, so two-point geostatistical assumptions are made regarding the reservoir properties and the initial ensemble is generated as random samples drawn from a multivariate Gaussian distribution defined by a variogram or covariance function. If historical data of nearby reservoirs or already induced fractures is available, such data can be used as a starting point for making assump-

9



tions regarding the mean and covariance function for generating the initial ensemble of realizations of the unknown reservoir properties. Such data can also be available, to an extent, from the minifrac test performed prior to the hydraulic fracturing operation. Using this a priori knowledge, an initial ensemble of R realizations is constructed by adding R different perturbations to the mean of the spatially varying parameter distribution. In the next step, called the assimilation step, the cross covariance between the predicted observations and the unknown parameters is updated by means of a filter gain using the available true measurements. In a hydraulic fracturing process, these available measurements are the estimates of fracture length and fracture width at the wellbore which are provided via the processed micro-seismic and downhole pressure data, respectively. Now, since the number of parameters to be estimated is huge compared to the available measurements, the inverse problem becomes ill-posed. One can increase the identifiability of this inverse problem by reducing the number of parameters to be estimated by transforming the original high-dimensional parameter space to a low-dimensional subspace. Also, since the spatially varying geological properties follow certain patterns, it is required to preserve the key spatial features observed in the high-dimensional parameter space. To achieve this, we apply POD to extract the most dominant geological features and provide a low-dimensional representation of the parameters, which is then incorporated in the EnKF parameter estimation routine. Below, we present a mathematical formulation for the steps involved in the POD-based EnKF estimation methodology. Let us consider a horizontal fracture geometry containing I grid points in the discretized system. Therefore, the representation of the spatially heterogeneous parameters (e.g. Young’s modulus profile) will be in the form of an I × 1 vector. One can generate R realizations of the spatially varying Young’s modulus profile as described above, and construct a representative ensemble of this data in the form of an I × R matrix. Now, we calculate the correlation matrix using the realizations and solve the integral eigenvalue problem to obtain the basis functions. The members of the

10


Page 10 of 43

Page 11 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ensemble can be approximately represented by a set of basis functions as follows: d

E ≈ ∑ µi ϕi

(10)

i=1

Here E ∈ RI×1 denotes the spatially varying Young’s modulus, which is an I-dimensional vector, ϕi ∈ RI denote the orthogonal basis functions, µi ∈ R are real-valued coefficients that define the appropriate linear combination of basis functions to approximate the true spatially varying Young’s modulus profile, and d is the dimension of the reduced parameter space obtained by examining the energy (in a mean square sense) of the singular values. The expression in Eq. (10) is a reducedorder representation of the parameter space where the basis functions, shared by all the ensemble members, contain the key spatial patterns of the unknown parameters. These basis functions induce a parametrization of the spatially varying Young’s modulus profile if we allow the coefficients µi to be free. Therefore, these coefficents will constitute the actual variables to be estimated via the EnKF, thereby effectively reducing the number of parameters to be estimated from I to d. Let the vector µ k ∈ Rd be the collection of all the free variables used in Eq. (10) at time tk . Now, at time step k, let us consider the forward (high-fidelity) model:

xk+1 = f(xk , uk , µ k )

(11)

where f is a nonlinear function given in Eqs. (1)-(9) of the state xk , the coefficients µ k and the input uk . Please note that we have used a subscript k for µ simply to denote that its value updates at every data assimilation step, but the model parameters do not change with time. Let yk ∈ R p be the predicted data (i.e., system outputs) at time tk . In this work, the measurements are the fracture length and the fracture width at the wellbore which are a subset of the state vector xk . The pdimensional vector of predicted outputs is therefore generated as yk = Axk , where A is a linear operator which contains zeroes and ones, and picks out the components of xk that are measured at time tk . However, since the states are a nonlinear function of the model parameters, the relation between the predicted outputs and the model parameters can be represented by a nonlinear operator, 11



Page 12 of 43

µ k , uk ). h, as yk = h(µ Since the measurements are a nonlinear function of the model parameters, we can consider an augmented state vector, Xk ∈ R(d+p)×1 , that includes the output measurements, yk along with the unknown model parameters, µ k : 







µ   µ  k  k       Xk =   =       µ k , uk ) yk h(µ

(12)

With this definition, the relation between the predicted output vector, yk , and the augmented state vector, Xk , is given by the following linear expression:  "



#   µk    I     µ k , uk ) h(µ

yk = HXk = 0

(13)

where H ∈ R p×(d+p) , 0 ∈ R p×d and I ∈ R p×p . In this work, the state variables, xk include the fracture width and proppant concentrations at various locations across the fracture and the fracture length, the model parameters include spatially varying Young’s modulus, although it is parameterized via POD and the coefficients, µ k , are considered as the actual variables to be estimated. Additionally, the observations, yk , include the fracture length and fracture width at the wellbore that are processed from the microseismic monitoring and downhole pressure analysis. The inputs to the system include the flow rate and the proppant concentration at the wellbore. Below is the sequence of steps involved in the estimation of spatially varying parameters using the POD-based EnKF scheme. ˆ = {E1 , E2 , . . . , ER } ∈ RI×R using 1. Generate R realizations of the unknown parameter fields E the mean and covariance functions available from the prior knowledge. 2. Parameterize the ensemble of realizations via POD, and compute the spatial basis functions 12


Page 13 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


µ 1k , µ 2k , . . . , µ Rk } ∈ Φ = {ϕ1 (z), ϕ2 (z), . . . , ϕd (z)} ∈ RI×d , leaving the coefficients, M = {µ Rd×R , to be free. 3. Given a vector of true measurements dk ∈ R p available at every sampling time tk , generate R artificial observations for data assimilation. j

j

j

dk = dk + εk , j = 1, . . . , R; εk ∈ N(0, CD )

(14)

which can be stored in the columns of Dk = {d1k , . . . , dRk } ∈ R p×R . Here ε is a random perturbation drawn from a normal distribution such that ε ∈ N(0, CD ), where CD ∈ R p×p is the covariance matrix of the measurement noise. 4. Forecast : At time tk , for every realization of the spatially varying parameters in the set f

1, f

R, f

M, the forecast of the ensemble of augmented states, Xk = {Xk , . . . , Xk } ∈ R(d+p)×R , is generated according to:  j, f Xk







j,a µ j, f    µ k−1  k    =  =      j, f    j,a µ k−1 , uk−1 ) yk h(µ

(15)

where the superscripts f and a denote the forecast and assimilation steps, respectively. 5. Assimilation: The covariance matrix at time tk , Ck ∈ R(d+p)×(d+p) , associated with the ensemble of augmented states is computed as: ¯f=1 X k R

Ck =

R

∑ Xkj, f

(16)

j=1

1 f ¯ f · 1R X f − X ¯ f · 1R T Xk − X k k k R−1

(17)

¯ f denotes the ensemble average of the augmented states and 1R ∈ R1×R is a row where X k vector whose elements are all equal to 1. The update to the predicted state and outputs is 13



Page 14 of 43

then computed as: Kk = Ck HT HCk HT + CD f

−1

(18)

f

Xak = Xk + Kk Dk − HXk

(19)

where Kk ∈ R(d+p)×R is the Kalman gain of the ensemble at time tk . 6. The sequence of forecast and assimilation steps is repeated until all the available measurements are assimilated. 7. The mean of the free coefficients in the final ensemble, µ¯ a ∈ Rd , along with the previously obtained basis functions, Φ , is used to compute the unknown spatially varying Young’s modulus profile, Etrue , as follows: d

Etrue ≈ ∑ µ¯ ia ϕi = Φ µ¯ a ,

µ¯ a =

i=1

1 R

R

∑ µ j,a

(20)

j=1

We refer the readers to Ref. 48 for rigorous details regarding the practical EnKF implementation.

Optimal pumping schedule design In this section, we focus on the design of a model-based feedback controller for the hydraulic fracturing process. It is presented as a two-step process: (1) estimation of the spatially varying Young’s modulus profile by using the process measurements via the POD-based EnKF methodology and (2) design of a model predictive controller by using reduced-order models computed using the local in time model order-reduction methodology developed in our previous work. 24 The schematic illustration of the overall methodology is presented in Fig. 2. The algorithm for the parameter estimation via POD-based EnKF is presented in the previous section. For the sake of brevity, the algorithm for generating reduced-order models for the purpose of controller design is not presented in this manuscript. However, the readers can refer to our previous work 24 for a detailed algorithm on the closed-loop control of the hydraulic fracturing process. 14


Page 15 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Within this context, we first perform numerical experiments for parameter estimation where the mean of the updated ensemble provides an accurate approximation for the true spatially varying Young’s modulus profile. In this work, we estimate only the spatially varying Young’s modulus and assume all other process model parameters are known. Following the parameter estimation, we solve the high-fidelity model of hydraulic fracturing with the updated spatially varying Young’s modulus profile to generate the data required for model order-reduction. The developed reducedorder models are then used to design the optimal pumping schedule required to achieve a uniform concentration profile across the fracture at the end of pumping. As mentioned earlier, the trapped proppant inside the fracture increases fluid conductivity by creating conduits through which oil and gas can be transported easily (because of high permeability of the proppant compared to the surrounding rock formation) to the wellbore where it can be extracted. Therefore, achieving uniform final proppant concentration to generate a uniform conductive fracture channel is vital in enhancing the productivity of a fractured well. There are several constraints associated with the hydraulic fracturing process that need to be accounted for in the controller design. We employ model predictive control schemes to handle these constraints in a systematic way in the design of the feedback control system. The following subsections provide a brief description regarding the reduced-order models and the model predictive control formulation employed in this work.

Reduced order models Due to the infinite-dimensional nature of the forward model, Eqs. (1)-(9), it cannot be directly used in the controller design. Alternatively, we seek computationally efficient reduced-order models generated by employing the LDMDc methodology developed in our previous work. 24 LDMDc is an equation-free model reduction framework that extracts the inherent system dynamics by distinguishing them from the applied external inputs. Moreover, the LDMDc framework results in linear discrete state-space models that can be directly used in the model predictive control formulation. The LDMDc is able to generate reduced-order models that provide accurate approximations to the high-fidelity solutions by effectively capturing the local dynamics in moving boundary problems

15



Page 16 of 43

such as hydraulic fracturing. By applying the LDMDc framework with the high-fidelity simulation data, we obtain the set of local in time reduced-order models of the following form: ˆ i xˆ (tk ) + B ˆ i u(tk ), xˆ (tk+1 ) = A

tk ∈ C i

yˆ (tk ) = Hˆx(tk ) C i = x j |(kx j − ci k22 ) ≤ (kx j − c` k22 ), ∀` 6= i , i = 1, . . . , m m [

(21)

C i = {x1 , . . . , xn }

i=1

ˆ i ∈ Rr×r and Bˆ i ∈ Rr×1 where xˆ ∈ Rr×1 denotes the vector of states in the reduced-order model, A denote the local in time process models corresponding to ith temporal subdomain, C i , and yˆ (tk ) denotes the vector of output variables. In the above equation x j , j = 1, . . . , n denote the state snapshots which are partitioned into m clusters whose centers are represented by ci , i = 1, . . . , m. We refer the readers to Ref. 24,43,49 for a comprehensive description of the domain partitioning problem. For the design of controllers, we assumed that the available real-time measurements include the proppant concentration C(tk ) at 6 different positions across the horizontal fracture along with the fracture width at the wellbore, W0 (tk ) and fracture length, L(tk ). However, in practice, measuring the proppant concentration inside the fracture is not feasible. In that case we can employ state estimators such as Kalman filter 6,50 to estimate the values of unmeasurable states.

Model predicitve control Once the reduced-order models are generated, we formulate the following MPC problem to design an optimal puming schedule required to achieve uniform concentration of the proppant across the

16


Page 17 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fracture at the end of pumping.

Minimize C0,k

s.t

ˆ f ) −Ctarget )T Qc (C(t ˆ f ) −Ctarget ) (C(t

(22a)

Reduced-order models, Eq. (21)

(22b)

xˆ (tk ) = x(tk )

(22c)

ˆ k + j∆) ≤ Cmax , Cmin ≤ C(t C0,k−1+m ≤ C0,k+m ≤ C0,k−1+m + 4, 2Q0 ∆ ∑ C0,k = M prop

∀ j = 0, . . . , 10 − k

(22d)

m = 1, . . . , 10 − k

(22e) (22f)

k

¯ f ) = Lopt , L(t

W¯ 0 (t f ) ≥ Wopt

(22g)

ˆ indicates the predicted state trajectory, Ctarget is the desired set-point value for the final where (·) proppant concentration, Qc is a positive definite matrix used to compute the weighted norm, t f ˆ is the vector of denotes the total treatment time, ∆ is the sampling time, tk is the current time, C proppant concentrations at 6 positions and C0,k is the inlet proppant concentration (i.e., manipulated input) corresponding to kth time interval i.e., t ∈ [tk ,tk+1 ). In the above optimization problem, the objective function, Eq. (22a), is a quadratic penalty that is the squared deviation of the proppant concentration from the set-point value at 6 different positions at the end of pumping. xˆ (t) denotes the predicted state trajectory computed using an approximate linear model, Eq. (22b), which is available to us by applying the LDMDc methodology as detailed in the previous section. Please note that since the developed reduced-order models are local in time, we have to use different models at different time instants. Therefore, the corresponding approximate model, Eq. (22b), is used to predict the state trajectory at a particular time interval ˆ i and B ˆ i main which the calculation is being performed. More specifically, at any time tk ∈ C i , A trices corresponding to the ith cluster (i.e., ith temporal subdomain) will be selected to predict the state trajectory. The initial conditions are given in Eq. (22c) which are obtained at each sampling time from the current full-state measurements. To account for the fact that the concentration of the

17



proppant inside the fracture should not exceed a certain limit at any time instant to avoid premature termination of the process, Eq. (22d) is included as a constraint. Eq. (22e) demands a monotonic increase in the input concentration profile with a maximum increment of 4 ppga/stage, where ppga is a concentration unit used in petroleum engineering that refers to 1 pound of proppant added to a gallon of water. A material constraint, Eq. (22f), is also introduced to ensure the required amount of proppant, M prop , is injected into the fracture. To make sure that the fracture achieves a desired geometry at the end of pumping, Eq. (22g) imposes a set of terminal constraints on the fracture geometry. Please note that the performance index in the above MPC formulation does not include a penalty term on the process inputs. Instead, this is implicitly addressed by introducing the material constraint, Eq. (22f).

Results Simulation description For the high-fidelity simulation, we generate a synthetic spatially varying Young’s modulus profile which is shown in Fig. 3 (black). This is taken to be the true parameter profile. Eqs. (1)-(9) constitute the high-fidelity model and they are solved using a higher-order discretization scheme to describe the fracture propagation, fluid flow, and proppant transport simultaneously. The fracture propagation is terminated at 135 m and the spatial domain is discretized with each grid point having a size 0.1 m, resulting in a total of I = 1351 points. The fluid injection rate is fixed at 0.03 m3 /s. The fracture width at the wellbore and the fracture length obtained from the simulation are considered as the true process measurements, giving p = 2. The values of the key process parameters used in the simulation are summarized in Table 1. To estimate the true Young’s modulus distribution generated above, we construct an initial ensemble of R = 100 parameter realizations. To achieve this, we consider a mean that is different from the true parameter profile to represent the initial distribution. In practice, as noted earlier, the first two statistical moments of the spatially varying Young’s modulus profile are known a priori 18


Page 18 of 43

Page 19 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Simulation parameters Property

Value

Fracture height, H - m

20

Leak-off coefficient, Cleak - m/s0.5

6.3 × 10−5

Proppant particle density, ρsd - kg/m3

2648

Pure fluid density, ρ f - kg/m3

1000

Fracture fluid viscosity, η - Pa · s

0.56

Poisson ratio of formation, ν

0.2

based on available historical data. The mean of the initial ensemble is shown in Fig. 3 (green). A covariance function is then used to generate 100 parameter realizations by adding 100 different perturbations to the initial ensemble mean. Fig. 3 (blue) shows two representative realizations selected from the initial ensemble that are at the extremes of the ensemble range. Each realization consists of I = 1351 points describing the spatial variability in the Young’s modulus profile. The high-fidelity model of the hydraulic fracturing process is solved with each realization that represents a spatially varying Young’s modulus profile to predict the fracture width at the wellbore and the fracture length. The predicted outputs are compared against the true measurements at every 25 s. The total duration of the process is 1000 s, which gives k = 40 time steps. Consequently, every parameter realization utilizes a set of 2 × 40 distinct measurements. In the numerical experiments, measurement errors are modeled as white noise (Gaussian random samples with zero mean) with the covariance CD , the value of which is taken as diag(0.005, 0.008), where the first term in the matrix corresponds to the error in the wellbore width and the second term to the fracture length, respectively. The units of the above matrix CD are meters and the above values have been determined based on the order of the wellbore width whose maximum value is aproximately 15 mm (the variance of the error is taken as 5 mm). During the forward simulations, the predicted outputs are assimilated against the true measurement values with added noise to update the ensemble at each time step. Therefore, we performed a total of 100 × 40 high-fidelity simulations.

19



Parameter estimation We first performed POD on the initial ensemble where 100 different spatially varying Young’s modulus profiles are sampled from the known first two statistical moments. The computed spatial basis functions are then employed to reduce the number of parameters to be estimated. Fig. 4 shows the cumulative energy content of the singular values associated with the ensemble. Based on this we chose the first d = 10 basis functions which correspond to 99% of the total energy captured µ ∈ R10×100 ) that describe the linear combination of by the singular values. The coefficients (µ these basis functions to reconstruct the parameter profiles are allowed to be free in the EnKF algorithm which leads to the parameterization of the spatially varying Young’s modulus profile. We emphasize that the number of parameters to be estimated in this problem was reduced to 10 from 1351, which alleviates the ill-posed nature of the original problem. Remark 2 Please note that the reduction in the number of parameters via POD depends on the covariance function chosen to sample the initial ensemble. This is because the rate of decay of singular values depends on the covariance function used. Once the parameterization step is completed, we implemented the EnKF algorithm to obtain the best estimate of the true spatially varying parameter profile. Since a major computational requirement is associated with the high-fidelity simulations used in the forecast step, the EnKF algoriothm was conveniently parallelized to speed up the computations. In our work, the highfidelity simulations and the EnKF programs are developed using Parallel Computing ToolboxTM R in MATLAB and then scaled up to computing clusters by running it on MATLAB Distributed

Computing ServerTM . Each assimilation step gives a set of new coefficient vectors that are used with the basis functions to regenerate the estimated spatially varying Young’s modulus profile. This set of new parameters is further used in the forecast of output variables at the next time step. The forecast and the assimilation steps are repeated until all the measurements are integrated into the high-fidelity simulation. Remark 3 We note that ill-conditioning might occur when computing the matrix inverse in the 20


Page 20 of 43

Page 21 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Kalman gain formulation, in Eq. (18). This potential singularity of the inverse computation requires the use of pseudo inverse via Truncated Singular Value Decomposition (TSVD). We would like to note that we have followed the practical implementation suggested by Ref. 48,51,52 in computing the pseudo inverse and using the singular vectors to complete the analysis step in the EnKF algorithm. The results of parameter estimation are shown in Fig. 5. Specifically, Fig. 5(a) shows an ensemble member, which represents a candidate Young’s modulus profile, after the final update step and Fig. 5(b) shows the final ensemble mean. In a typical EnKF implementation, the mean of the final ensemble is taken as the approximation for the unknown parameter profile. Therefore, it is evident from the figure that the POD-based EnKF method provided a very good estimate for the true spatially varying Young’s modulus profile after assimilating all the available measurement data. The estimated parameter profiles in Fig. 5 retain the major spatial trend of the true Young’s modulus profile. This is attributed to the basis functions which capture the dominant spatial patterns in the underlying data of the prior distribution. Fig 6 shows the evolution of error covariances during the data assimilation process. One would expect the covariances to approach to zero as more data is being assimilated. This can be observed in Fig 6 where both variances (Cii ) and covariances (Ci j ) decline to zero with time. Additionally, this also shows that the ensemble estimates converge to the true values at the end of data assimilation. Although, the figures show a few elements taken from the covariance matrix to describe their time evolution, a similar behavior is observed throughout the covariance matrix. Please note that the covariances have been normalized prior to the plots. One way to measure the performance of the EnKF method is the reduction in the final ensemble spread which is evident by comparing Figs. 5(a) and 5(b). From the figure, it can be seen that the difference between the final ensemble member and the mean is small. Note that for illustration purposes, only one representative member is selected randomly from the final ensemble. However, this reduction in the spread is observed with respect to all the final ensemble members. Furthermore, by comparing Figs. 3 and 5, it can be observed that the uncertainty in the ensemble members is also reduced between the prior and posterior distributions. In order to highlight this point more 21



clearly, we show the ratio of the standard deviation between the final and initial ensemble in Fig. 7. Please note that a value close to 0 indicates a large reduction in uncertainty whereas a value close to 1 indicates a small reduction. As can be seen from the figure, a significant reduction in uncertainty is achieved throughout the fracture. Fig. 8 shows the mean square error (MSE) plot between the true parameter values and the mean of the ensemble posterior distribution at each time step. Note that the MSE values have been normalized with respect to the maximum. As can be seen from the figure, a significant reduction in the error is observed between the initial and final update steps. As more data are integrated into the high-fidelity model, the EnKF estimate converges to the original value, resulting in a decrease in the MSE . In addition to finding a good approximation between the true and estimated parameter values, one of the primary objectives of parameter estimation of dynamic models is to accurately predict the process dynamics. This is an important feature in order to design model-based control systems or perform dynamic optimzations of the process. The comparison between the process dynamics predicted by the updated model, based on the estimated parameters, and the true values is shown in Fig. 9. The predictive capabilities of the updated model is evident from the figure which shows a good match between the true and predicted values of the fracture width at the wellbore and the fracture length through Fig. 9(a) and 9(b), respectively. Remark 4 Please note that although the EnKF handles the nonlinear dynamics correctly during the forecast step, it sometimes fails in the analysis step. 53–56 More specifically, for a highly nonlinear problem, the current state estimates obtained from EnKF and the current state values computed by re-solving the forward model with the updated parameters may not be identical, implying the production of unrealistic parameter estimates. However, we emphasize that in this manuscript we only estimated the model parameters via the EnKF and the predicted outputs are obtained from the forward model with the updated model parameters. Remark 5 We would like to note that, when we included the states into the augmented state vector of Eq. (12), the above described inconsistency was observed and the EnKF gave unrealistic values for the current state estimates, the results of which have not been reported in this manuscript. 22


Page 22 of 43

Page 23 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Closed-loop control With the updated spatially varying Young’s modulus profile, open-loop simulations of the highfidelity model are performed using a set of training inputs to generate the data required for developing the reduced-order models. The LDMDc methodology determined the optimal number of clusters as m = 100 and the corresponding reduced-order models are generated within each subtemporal domain. Fig. 10 shows the comparison between the actual state values and their corresponding approximate solutions of the fracture width at the wellbore, fracture length, and proppant concentration profiles at two different positions. It can be observed that the LDMDc technique provides accurate approximations to the true process state values obtained from the high-fidelity model. Now that we have accurate reduced-order models, a model-based feedback controller is designed to achieve a uniform final proppant concentration profile across the fracture along with a desired final fracture length and width to enhance the production of the fractured well. The values used for various process parameters are the same as the ones given in Table 1. In the closed-loop simulation, the sampling time ∆ and the total time of operation t f were chosen to be 100 s and 1220 s, respectively, during which a total of M prop = 48, 000 kg of proppant has been injected into the fracture. The sampling time was chosen to be around 100 s to accommodate the practical time required for processing the raw microseismic data and downhole pressure analysis to obtain the wellbore fracture width and fracture length measurements, respectively. During initial stages of the hydraulic fracturing process, a high-pressure fluid (called pad) consisting mostly of water is pumped to break the rock and propagate fractures in the formation at perforated sites. The duration of this pad stage, t p , has been fixed at 220 s to reach the desired fracture length of Lopt = 135 m without tip-screen out (premature termination of the process). The feedback control system is initialized after the injection of pad (i.e., tk ≥ t p ). The proppant pumping schedule is partitioned into 10 substages (i.e, k = 1, . . . , 10) and the duration of each substage is given by the sampling time, ∆. We assumed that at the beginning of each substage, the full-state measurements of x(tk ) are available (x = [W0 , L, C]), which are used to predict the future states via the developed local in time 23



reduced-order models. The predicted state trajectories are then used in the optimization problem to compute the control inputs, and the corresponding process behavior that minimizes the squared deviation from its prespecified target value. The uniform target proppant concentration desired at the end of pumping is Ctarget = 10 ppga. Please refer to Ref. 6 for details on the calculation of this set-point value. The positive definite matrix, Qc , containing the weights on the penalty function of the MPC optimization problem is considered to be an identity matrix. However, one can use appropriate weights depending on the process nature. The first step of the input profile, C0,k , obtained by solving the optimization problem over a prediction horizon length of N p was applied to the high-fidelity model in a sample-and-hold fashion, and this procedure was repeated at every sampling time until the end of pumping. In this application, we used a shrinking horizon which gives N p = t f − tk . We solve the optimization problem to compute control inputs while regulating the output value to be close to the set-point at the end of the final time step, which is what we desire in this specific application. The results are presented in Fig. 11(a) which shows the final concentration profile of the closedloop system under MPC. Clearly, the designed controller drives the final proppant concentration (red) towards the target value (black) while propagating the fracture to the desired length (arrow) specified by the geometry constraints. For the sake of illustrating the effect of not characterizing the spatially varying Young’s modulus profile prior to the design of feedback control systems, we have also designed a closed-loop system under MPC based on a constant Young’s modulus profile. Here, the value of the Young’s modulus is fixed at E = 50 × 102 MPa across the fracture. The closed-loop performance under the MPC designed based on the constant Young’s modulus is shown in Fig. 11(b). As can be seen, when the rock propertiers are not characterized accurately, it can result in a poor controller performance, i.e, non-uniform final concentration profiles, thereby diminishing the production capabilities of the fracture well.

24


Page 24 of 43

Page 25 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


CONCLUSIONS A combined POD based EnKF method for the estimation of spatially varying Young’s modulus profile for the closed-loop control of a hydraulic fracturing process to achieve uniform proppant concentration at the end of pumping is presented. Parametrizing via POD resulted in a significant increase in the identifiability of the original parameter estimation problem, i.e., the number of variables to be estimated in the data assimilation framework is reduced from 1351 to 10. The method gives good results after all the available data has been integrated, and it is demonstrated that a significantly better match is obtained for the Young’s modulus estimate relative to the synthetic true Young’s modulus profile that is used for testing purposes. Additionally, we observed significant predictive capabilities from the updated model as evident from the close match to both the fracture length and the fracture width at the wellbore. The designed closed-loop controller based on the updated parameter profile, is able to drive the proppant concentration profile at each spatial location to a uniform target value at the end of pumping. Comparison of closed-loop simulation results shows the significance of accurately characterizing the spatial variability in the reservoir parameters.

25



References 1. Economides, M. J.; Martin, T. Modern fracturing: Enhancing natural gas production; Energy Tribune Publishing: Houston, 2007. 2. Economides, M. J.; Nolte, K. G. Reservoir stimulation; Wiley: Chichester, 2000. 3. Gu, H.; Desroches, J. New pump schedule generator for hydraulic fracturing treatment design. SPE Latin American and Caribbean Petroleum Engineering Conference, Port-of-Spain, Trinidad and Tobago, April 27-30, 2003. 4. Dontsov, E. V.; Pierce, A. P. A new technique for proppant schedule design. Hydraulic Fracturing J. 2014, 1, 1–8. 5. Gu, Q.; Hoo, K. A. Model-based closed-loop control of the hydraulic fracturing process. Ind. Eng. Chem. Res. 2015, 54, 1585–1594. 6. Siddhamshetty, P.; Yang, S.; Kwon, J. S. Modeling of hydraulic fracturing and designing of online pumping schedules to achieve uniform proppant concentration. Comput. Chem. Eng. In Press 2017, DOI: https://doi.org/10.1016/j.compchemeng.2017.10.032. 7. Siddhamshetty, P.; Liu, S.; Valko, P.; Kwon, J. S. Feedback control of proppant bank heights during hydraulic fracturing for enhanced productivity in shale formations. AIChE J. In press 2017, DOI: 10.1002/aic.16031. 8. Daniels, J. L.; Waters, G. A.; Le Calvez, J. H.; Bentley, D.; Lassek, J. T. Contacting more of the barnett shale through an integration of real-time microseismic monitoring, petrophysics, and hydraulic fracture design. SPE Annual Technical Conference and Exhibition, Anaheim, California, November 11-14, 2007. 9. King, G. E. Thirty years of gas shale fracturing: What have we learned? SPE Annual Technical Conference and Exhibition, Florence, Italy, September 19-22, 2010.

26


Page 26 of 43

Page 27 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10. Evensen, G. Sequential data assimilation with a non-linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res. 1994, 99, 10143–10162. 11. Evensen, G. Advanced data assimilation for strongly nonlinear dynamics. Mon. Weather Rev. 1997, 125, 1342–1354. 12. Park, J. H.; Kaneko, A. Assimilation of coastal acoustic tomography data into a barotropic ocean model. Geophys. Res. Lett. 2000, 27, 3373–3376. 13. Haugen, V. E.; Evensen, G. Assimilation of SLA and SST data into an OGCM for the Indian ocean. Ocean Dynamics 2002, 52, 133–151. 14. Mitchell, H. L.; Houtekamer, P. L.; Pellerin, G. Ensemble size, and model-error representation in an ensemble kalman filter. Mon. Weather Rev. 2002, 130, 2791–2808. 15. Warpinski, N. R.; Engler, B. P.; Young, C. J.; Peterson, R.; Branagan, P. T.; Fix, J. E.; James, E. Microseismic mapping of hydraulic fractures using multi-level wireline receivers. SPE Annual Technical Conference and Exhibition, Dallas, Texas, October 22-25, 1995. 16. Rutledge, J. T.; Phillips, W. S. Hydraulic stimulation of natural fractures as revealed by induced microearthquakes, Carthage Cotton Valley gas field, east Texas. Geophysics 2003, 68, 441–452. 17. Maxwell, S. C.; Rutledge, J.; Jones, R.; Fehler, M. Petroleum reservoir characterization using downhole microseismic monitoring. Geophysics. 2010, 75, 129–137. 18. Holmes, P.; Lumley, J. L.; Berkooz, G. Turbulence, Coherent Structures, Dynamical Systems and Symmetry; Cambridge University Press: New York, 1996. 19. Berkooz, G.; Holmes, P.; Lumley, J. The proper orthogonal decomposition in the analysis of turbulent flows. Ann. Rev. Fluid. Mech. 1993, 25, 539–575.

27



20. Sarma, P.; Durlofsky, L. J.; Aziz, K.; Chen, W. H. Efficient real-time reservoir management using adjoint-based optimal control and model updating. Computational Geosciences 2006, 10, 3–36. 21. Arfa, S.; Gildin, E. Permeability parametrization using higher order singular value decomposition. Proceedings of the 12th International Conference on Machine Learning and Applications, Miami, Florida, December 4-7, 2013; pp 188-193. 22. Arfa, S.; Gildin, E. Tensor based geology preserving reservoir parameterization with Higher Order Singular Value Decomposition (HOSVD). Comput. Geosci. 2016, 94, 110–120. 23. Jafarpour, B. Sparsity-promoting solution of subsurface flow model calibration inverse problems. Adv. Hydrogeol. 2013, 73–94. 24. Narasingam, A.; Kwon, J. S. Development of local dynamic mode decomposition with control: Application to model predictive control of hydraulic fracturing. Comput. Chem. Eng. 2017, 106, 501–511. 25. Perkins, T. K.; Kern, L. R. Widths of hydraulic fractures. J. Pet. Technol. 1961, 13, 937–949. 26. Gu, Q.; Hoo, K. A. Evaluating the Performance of a Fracturing Treatment Design. Ind. Eng. Chem. Res. 2014, 53, 10491–10503. 27. Yang, S.; Siddhamshetty, P.; Kwon, J. S. Optimal pumping schedule design to achieve a uniform proppant concentration level in hydraulic fracturing. Comput. Chem. Eng. 2017, 101, 138–147. 28. Howard, G. C.; Fast, C. R. Optimum fluid characteristics for fracture extension. Drill. prod. pract. 1957, 24, 261–270. 29. Adachi, J.; Siebrits, E.; Pierce, A.; Desroches, J. Computer simulation of hydraulic fractures. Int. J. Rock Mech. Min. Sci. 2007, 44, 739–757.

28


Page 28 of 43

Page 29 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


30. Daneshy, A. Numerical solution of sand transport in hydraulic fracturing. J. Pet. Technol. 1978, 30, 132–140. 31. Barree, R.; Conway, M. Experimental and numerical modeling of convective proppant transport. J. Pet. Technol. 1995, 47, 216–222. 32. Novotny, E. J. Proppant transport. SPE Annual Fall Technical Conference and Exhibition, Denver, Colorado, October 9-12, 1977. 33. Liu, E.; Tod, S. R.; Li, X. Y. The effects of stress and pore fluid pressure on seismic anisotropy in cracked rocks. Canadian SEG Recorder 2002, 27, 92–98. 34. Tod, S. R. The effects of stress and fluid pressure on the anisotropy of interconnected cracks. Geophys. J. Int. 2002, 149, 149–156. 35. Maultzsch, S.; Chapman, M.; Liu, E.; Li, X. Y. Modelling frequency dependent seismic anisotropy in fluid-saturated rock with aligned fractures: implicaition of fracture size estimation from anisotropic measurements. Geophys. Prosp. 2003, 51, 381–392. 36. Liu, E.; Queen, J. H.; Li, X. Y.; Chapman, M.; Maultzsch, S.; Lynn, H. B.; Chesnokov, E. M. Observation and analysis of frequency-dependent seismic anisotropy from a multi-component VSP. J. Appl. Geophys. 2003, 54, 319–333. 37. Zimmerman, R. W.; Bodvarsson, G. S. Hydraulic conductivity of rock fractures. Transp. Porous Media 1996, 23, 1–30. 38. Oron, A. P.; Berkowitz, B. Flow in rock fractures: the local cubic law assumption reexamined. Water Resour. Res. 1998, 34, 2811–2825. 39. Renshaw, C. E. On the relationship between mechanical and hydraulic apertures in roughwalled fractures. J. Geophys. Res. 1995, 100, 24629–24636.

29



40. Sisavath, S.; Al-Yaaruby, A.; Pain, C. C.; Zimmerman, R. W. A simple model for deviations from the cubic law for a fracture undergoing dilation or closure. Pure Appl. Geophys. 2003, 160, 1009–1022. 41. Méheust, Y.; Schmittbuhl, J. Scale effects related to flow in rough fractures. Pure Appl. Geophys. 2003, 160, 1023–1050. 42. Zimmerman, R. W.; Main, I. In Mechanics of Fluid-Saturated Rocks; Gueguen, Y., Bouteca, M., Eds.; Elsevier Academic Press: London, 2004; Chapter Hydromechanical behaviour of fractured rocks, pp 361–419. 43. Narasingam, A.; Siddhamshetty, P.; Kwon, J. S. Temporal clustering for order reduction of nonlinear parabolic PDE systems with time-dependent spatial domains: Application to a hydraulic fracturing process. AIChE J. 2017, 63, 3818–3831. 44. Burgers, G.; van Leeuwen, P. J.; Evensen, G. Analyisis scheme in the ensemble Kalman Filter. Mon. Weather Rev. 1998, 126, 1719–1724. 45. Li, G.; Han, M.; Banerjee, R.; Reynolds, A. C. Integration of well test pressure data into heterogeneous geological reservoir models. SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, October 4-7, 2009. 46. Gao, G.; Zafari, M.; Reynolds, A. C. Quantifying the uncertainty for the PUNQ-S3 problem in a Bayesian setting with RML and EnKF. SPE J. 2006, 11, 506–515. 47. Liu, E. Effects of fracture aperture and roughness on hydraulic and mechanical properties of rocks: implication of seismic characterization of fractured reservoirs. J. Geophys. Eng. 2005, 2, 38–47. 48. Evensen, G. The Ensemble Kalman Filter: theoretical formulation and practical implementation. Ocean Dynamics 2003, 53, 343–367.

30


Page 30 of 43

Page 31 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


49. Tan, M. P.; Broach, J. R.; Floudas, C. A. A novel clustering approach and prediction of optimal number of clusters: global optimum search with enhanced positioning. J. Glob. Optim. 2007, 39, 323–346. 50. Kalman, R. E. A new approach to linear filtering and prediction problems. J. Basic Eng. 1960, 82, 35–45. 51. Skjervheim, J. A.; Evensen, G.; Aanonsen, S. I.; Ruud, B. O.; Johansen, T. A. Incorporating 4D seismic data in reservoir simulation models using ensemble Kalman filter. SPE J. 2007, 12, 282–292. 52. Rafiee, J.; Reynolds, A. C. Theoretical and efficient practical procedures for the generation of inflation factors for ES-MDA. Inverse Problems 2017, 33, 115003. 53. Thulin, K.; Li, G.; Aanonsen, S. I.; Reynolds, A. C. Estimation of initial fluid contacts by assimilation of production data With EnKF. SPE Annual Technical Conference and Exhibition, Anaheim, California, November 11-14, 2007. 54. Gu, Y.; Oliver, D. S. An iterative ensemble Kalman filter for multiphase fluid flow data assimilation. SPE J. 2007, 12, 438–446. 55. Zafari, M.; Reynolds, A. C. Assessing the Uncertainty in Reservoir Description and Performance Predictions With the Ensemble Kalman Filter. SPE J. 2007, 12, 382–391. 56. Wang, Y.; Li, G.; Reynolds, A. C. Estimation of depths of fluid contacts by history matching using iterative ensemble-Kalman smoothers. SPE J. 2010, 15, 509–525.

31



(a) Spatio-temporal evolution of the fracture width.

(b) Spatio-temporal evolution of the proppant concentration.

Figure 1: Profiles showing the spatio-temporal evolution of the states in a hydraulic fracturing process. It is important to note that the zero values in the above figures indicate the fracture has not yet propagated to the location. 32 ACS Paragon Plus Environment

Page 32 of 43

Page 33 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: Overall flow integrating the POD-based EnKF parameter estimation scheme with LDMDc-based MPC.

33



120 100 E(z), x 102 MPa

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

80 60 40 20 0

0

50

100 z, m

Figure 3: The true Young’s modulus profile (black) along with two representative realizations (blue) and the mean (green) of the initial ensemble.

34


Page 34 of 43

Page 35 of 43

1 % Energy

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.995

0.99 0

20

40 60 80 Number of basis functions

100

Figure 4: Cumulative energy captured by the basis functions as represented by their normalized singular values.

35



E(z), x 102 MPa

80 60 40 20 0 0

50

100

150

z, m (a) Final ensemble mean (red) and the true parameter profile (black)

80 E(z), x 102 MPa

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 43

60 40 20 0 0

50

100

150

z, m (b) A realization from the final ensemble (blue) and the true profile (black).

Figure 5: A comparison between the true and the final estimated Young’s modulus profile through POD-based EnKF. 36


Page 37 of 43

1 C11 C22

0.8

C33 C44 C55

Cii(t)

0.6

0.4

0.2

0 0

10

20 time step, k

30

40

(a) Evolution of the error variances, Cii , during data assimilation via the EnKF.

1 C13 C26

0.8

C37 C510 C78

0.6 Cij(t)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.4

0.2

0 0

10

20 time step, k

30

40

(b) Evolution of the error covariances, Ci j , during data assimilation via the EnKF.

Figure 6: Profiles showing the evolution of the error covariances, Ck , during the data assimilation steps. 37



1 STD ratio

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 43

0.5

0

0

20

40

60 80 z, m

100

120

Figure 7: Posterior to prior standard deviation ratio of the spatially varying Young’s modulus profile.

38


Page 39 of 43

1

0.8

0.6 MSE

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.4

0.2

0

0

10

20 time step, k

30

Figure 8: The profile of MSE with time.

39


40


W0(t), m

0.015 0.01 0.005 0

predicted true

0

200

400 600 time, s

800

1000

(a) Fracture width at the wellbore, W0 (t), predicted by the updated model (red) and its true values (black)

150 L(t), m

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 43

predicted true

100 50 0

0

200

400 600 time, s

800

1000

(b) Fracture length, L(t), predicted by the updated model (red) and its true values (black)

Figure 9: Profiles of the predicted system outputs and their true measurements.

40


Page 41 of 43

140

0.02 high-fidelity LDMDc

100

L(t), m

W0(t), m

high-fidelity LDMDc

120

0.015

0.01

80 60 40

0.005

20 0

0 0

500

1000

1500

0

500

(a) W0 (t) at the wellbore

1500

(b) L(t) of the fracture

10

10 high-fidelity LDMDc

high-fidelity LDMDc

8

C(t), ppga

8 6 4

6 4 2

2 0

1000

time, s

time, s

C(t), ppga

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

500

1000

0

1500

0

time, s

500

1000

1500

time, s

(c) C(t) at z = 18 m

(d) C(t) at z = 90 m

Figure 10: Comparison of the approximate solutions computed using local in time reduced-order models with the true measurements of the fracture width at the wellbore, fracture length and proppant concentration at two different positions.

41



12

C(tf), ppga

10 8 6 4 2 0

MPC - E(z) set-point

0

L = 135 m

50

100

z, m (a) MPC based on the updated spatially varying Young’s modulus E(z).

12 10

C(tf), ppga

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 43

8 6 4 2 0

MPC - E set-point

0

50

z, m

100

(b) MPC based on the constant Young’s modulus assumption.

Figure 11: Final proppant concentration profiles of the closed-loop system under MPC of hydraulic fracturing. 42


Page 43 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12: For Table of Contents Only

43


Handling spatial heterogeneity in reservoir parameters using POD

Recommend Documents