# «OPTIMIZATION AND MONITORING OF GEOLOGICAL CARBON STORAGE OPERATIONS A DISSERTATION SUBMITTED TO THE PROGRAM IN EARTH, ENERGY AND ENVIRONMENTAL ...»

## 4.3. DATA ASSIMILATION AND SENSOR PLACEMENT

4.3 Data assimilation and sensor placement We now describe our approach for performing data assimilation, which could be applied multiple times during closed-loop operation, as well as our procedure for placing monitoring wells or sensors. In our treatment, the sensors are placed after the injection wells are drilled but prior to any other data collection, so determining their locations is essentially an a-priori decision.

4.3.1 Data assimilation procedure A general approach for data assimilation involves stochastic inversion of observation data to ﬁnd a posterior probability density function (PDF) of the uncertain model parameters. Such a method is often framed within a Bayesian setting, and in practice entails deﬁning a prior distribution (e.g., in terms of geostatistical parameters) and then constructing samples from the posterior distribution by applying an inversion procedure. With this approach, one would seek to minimize the misﬁt function S

**given by Tarantola (2005):**

The ﬁrst term on the right hand side of Equation 4.4 quantiﬁes the mismatch between observed and simulated data, and the second (model mismatch) term provides a regularization that forces the history matched geological model to resemble the prior model. The speciﬁc terms are deﬁned as follows: m is a vector of unknown geological ¯ parameters (e.g., porosity and permeability in every grid block) with prior mean m and dimension of the order of the number of grid blocks in the model (dim(m) = O(Nb )), dobs ∈ RNd is a vector of observed data, where Nd is the number of observed data, f : RNb → RNd is the forward model, which maps input model parameters to output model data via ﬂow simulation, CD is the data covariance matrix, and CM is the covariance matrix deﬁning the spatial correlation structure of the prior model.

The representation in Equation 4.4 is strictly applicable only for Gaussian models (as are considered here). For systems characterized by multipoint spatial statistics, an alternate formulation may be required. There are a number of speciﬁc approaches that have been applied to perform data assimilation within the general context of

## 62 CHAPTER 4. CLOSED-LOOP AQUIFER MANAGEMENT

Equation 4.4, as discussed by Oliver et al. (2008).In this work, the unknown model parameters m refer to the porosity (φ) for each grid block of the model. We relate the x-component of grid block permeability (kx ) ¯ to porosity using Equation 2.12, where φ = 0.2, σφ = 0.025, a = 4.5 and b = 2 are speciﬁed parameters. As discussed in Chapter 2, the other permeability components are speciﬁed as ky = kx and kz = 0.1kx.

Our approach aims at ﬁnding a single (maximum likelihood) estimate of the geological model rather than generating samples from the posterior PDF. Finding the model m that minimizes an appropriately deﬁned misﬁt function can still be challenging due to the large number of unknown porosity variables in a typical simulation model. The number of optimization variables can be reduced considerably, however, through use of the Karhunen-Lo`ve (K-L) expansion (Loeve, 1977), as detailed later e in Section 4.3.2. Denoting the K-L parameters as ξ ∈ Rℓ, with ℓ Nb, we now express m as m(ξ). The use of this treatment allows us to remove the model mismatch term from the misﬁt function, as the K-L expansion incorporates the prior spatial correlation structure into the representation.

We also apply data weighting parameters to enhance the ability of the history matched model to predict CO2 plume location, which is the quantity of most interest here. Introducing a general weighting matrix W in place of C−1 in Equation 4.4, D as well as the K-L representation for the geological model m, the data assimilation

**problem becomes:**

ˆ min S = [f (m(ξ)) − dobs ]T W[f (m(ξ)) − dobs ], (4.5) ξ ˆ where S is the new misﬁt function to be minimized and f and dobs are again the modeled and observed data. The observed data here include time-dependent pressure data at observation wells, time-dependent bottomhole pressure (BHP) data at CO2 injectors, and the CO2 breakthrough time (BTT) at observation wells. We will discuss how we determine an optimal W in Section 4.3.4.

Because the K-L expansion provides a diﬀerentiable representation of the geological model m, the data assimilation problem in Equation 4.5 is suited to gradient based and direct search optimization methods that quickly converge to local minima.

## 4.3. DATA ASSIMILATION AND SENSOR PLACEMENT

See, e.g., Sarma et al. (2006) for the application of the K-L representation with gradient based optimization. Since we are using a commercial simulator and do not have access to gradients, we apply a local search procedure – Generalized Pattern Search (GPS, described previously in Section 2.2.3). We also investigated the use of PSO for the data assimilation problem, but we found GPS to consistently provide better matches than PSO for a given computational eﬀort. In addition, the more systematic local search performed by GPS (as opposed to the stochastic search using PSO) may be beneﬁcial in closed-loop applications. Speciﬁcally, GPS can use the history matched model from the previous update step as the initial guess, and the new model will then be in some sense ‘close’ to the previous model.4.3.2 Karhunen-Lo`ve representation of geological models e In this work we represent geological models using a discrete K-L expansion. This approach is also referred to as principal component analysis (PCA) or proper orthogonal decomposition (POD). Our K-L representation deﬁnes porosity φ, from which permeability kx is determined by the relationship in Equation 2.12. We consider models in which φ is normally distributed, which means kx is log-normally distributed. K-L representations are well suited for Gaussian models of this type.

**Using the K-L expansion for the heterogeneous porosity ﬁeld, we can write:**

** m(ξ) = Φξ + m, (4.6)**

where m ∈ RNb is the porosity ﬁeld (recall that Nb is the number of grid blocks in the model), Φ ∈ RNb ×ℓ is the basis matrix (Φ is a tall and skinny matrix, since ℓ Nb ), ξ ∈ Rℓ is the vector of geological parameters in the low-dimensional subspace, and m is the mean porosity ﬁeld. The K-L representation of the geological model has two main advantages: (1) realizations are deﬁned by a small set of parameters ℓ, which reduces the computation required for data assimilation, and (2) the resulting realizations of m(ξ) honor hard data such as porosity in well blocks and the prescribed correlation structure.

The basis matrix Φ is constructed by performing a singular value decomposition (SVD) of a set of prior realizations of m. Our speciﬁc approach, described in Sarma

## 64 CHAPTER 4. CLOSED-LOOP AQUIFER MANAGEMENT

et al. (2006), Sarma et al. (2008) and He et al. (2013), is as follows. Given a set of geostatistical (e.g., variogram) parameters and some amount of hard data, we apply SGeMS (Remy et al., 2009) to generate a set of Nr realizations of m. Here Nr is**typically O(102 − 103 ). The data matrix X ∈ RNb ×Nr is then given by:**

** X = [(m1 − m) (m2 − m)... (mNr − m)], (4.7)**

Nr where mi, i = 1... Nr, designates a particular realization and m = (1/Nr ) mi.

i=1 The columns of the basis matrix Φ are given by the left singular vectors of the SVD √ of X/ Nr, scaled by the corresponding principal values. There are a maximum of Nr nonzero singular values (and thus a maximum of Nr columns in Φ), though typically only ℓ columns are retained, where ℓ Nb. An energy criterion, with energy quantiﬁed in terms of the singular values associated with Φ, is often used to determine an appropriate value for ℓ (Sarma et al., 2006). In the models considered here, we have Nb = 5000 and we take Nr = 2000 and ℓ = 100. The construction of the K-L basis matrix represents a small computational overhead, as the SGeMS simulations are very fast and the SVD computation is not time consuming for Nr ∼ O(103 ).

New realizations that honor the prior spatial correlation structure and hard data can be constructed by generating standard normal ξ and then applying Equation 4.6.

Within the context of data assimilation, the geological model is represented using Equation 4.6 along with Equation 2.12. This means that, rather than having to determine Nb porosity values, we now must ﬁnd only the ℓ parameters appearing in the K-L expansion. This simpliﬁes considerably the history matching problem.

4.3.4 Optimal sensor placement and data weighting The data assimilation step entails modiﬁcation of the geological model in order to achieve agreement between ‘actual data’ and simulated results. The quality of the resulting geological models, in terms of their ability to predict quantities of interest, depends however on the speciﬁc data that are used in the history matching. In this section, we describe our procedure for determining the locations for sensors such that they are maximally eﬀective (in an optimization sense) in providing predictive models. We additionally determine optimal data weight parameters within the weighting matrix W, again with the intent of generating predictive models for which expected error is minimized.

The problem of placing sensors optimally is often referred to as Optimal Experimental Design (OED). Two types of methods are commonly applied to solve OED problems. The ﬁrst approach uses information theory, attempting to maximize the ‘information’ or ‘entropy’ in the expected data. The information in a dataset can be characterized in a number of diﬀerent ways (Haber et al., 2008). One of the most popular information measures is the Shannon Entropy, which was used by Magnant (2011) to determine the optimal placement of downhole gauges for cross-well seismic analysis in CO2 storage operations. The second approach for solving OED problems is to directly optimize the sensor placements for some measure of success. This method is often referred to as Bayesian experimental design because it typically entails characterizing a posterior PDF or maximum likelihood estimate for some measure of the

## 66 CHAPTER 4. CLOSED-LOOP AQUIFER MANAGEMENT

** Figure 4.2: Flowchart illustrating history matching based on GPS minimization and the K-L representation of geology **

**4.**

3. DATA ASSIMILATION AND SENSOR PLACEMENT

eﬀectiveness of the experimental design parameters.In this work, we use a simpliﬁed Bayesian experimental design algorithm (Chaloner & Verdinelli, 1995) to optimize sensor placements y and data weight parameters w (deﬁned below), which appear in W, to maximize the eﬃcacy of history matching. Sensor placements and data weight parameters are determined with the goal of minimizing an estimate of the error in the top layer CO2 mobility (as deﬁned in Equation 4.2) between the ‘true’ and history matched models. Since in practice we do not know the ‘true’ geological model, we minimize the expected value of mobility error by considering many ‘potentially true’ models. These potentially true models are consistent with the prior model and honor any data known at this point in the operation. Here, we assume that the porosity/permeability at the injection wells is known and we condition geologic realizations to that hard data.

Potentially true models are denoted mj, j ∈ 1... Ns. To each mj, we associate a history matched model, which we designate m∗. These m∗ can be viewed as the j j maximum likelihood estimate for mj following the application of the data assimilation procedure. The determination of optimal sensor placements and data weights thus amounts to minimizing the expected error (in top layer CO2 mobility) between potentially true models (mj ) and their corresponding history matched representations (m∗ ).

**This optimization problem can be formally stated as:**

j

ˆ with λi as deﬁned in Equation 4.2 and T = 500 years.

In order to evaluate F for any particular y and w, we need to have history matched models m∗ for each mj. This means that, for a single evaluation of F (for some y, w), j we must ﬁrst perform Ns history matches, and for each history match we then compute ∆[mj, m∗ (y, w)]. Becuase each history match may require thousands of simulations, j this direct approach would lead to a prohibitively expensive optimization. Therefore, we apply an alternative treatment.

Our approach is as follows. We generate a set of Ns geologic realizations (based on the prior), and then perform a ﬂow simulation for each of these models. In these simulations, injection well placements and rates are as determined from the a-priori optimization described earlier. The states (pressure and saturation) in every grid block at every time step in each of the Ns models are saved. Then, for each mj, rather than determine m∗ through detailed history matching, we assign m∗ to be the j j ˆ ‘closest’ of the Ns previously simulated models in terms of the data misﬁt S.

Our ‘approximate’ history matching procedure is to ﬁnd the model i ∈ {1,..., Ns }, i = j, such that the data mismatch with mj is minimized. This optimization problem

**is written as:**

ˆ

**where S is the data misﬁt function between models mi and mj :**

ˆ S(y, w, mj, mi ) = [f (mi, y) − dobs (mj, y)]T W(w)[f (mi, y) − dobs (mj, y)]. (4.11) Here f (mi, y) is the model data from realization mi with speciﬁed sensor positions y, and dobs is the observation data from realization mj with the same sensor positions as in the evaluation of f.