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

The optimization procedures applied in this work are noninvasive with respect to the ﬂow simulator and do not require gradient formulations. When computational resources are limited, we will mainly apply the Hooke-Jeeves Direct Search (HJDS) method (Hooke & Jeeves, 1961). HJDS has previously been shown to be eﬀective in optimizing subsurface ﬂow systems in both hydrological (Matott et al., 2011) and

## 18 CHAPTER 2. OPTIMIZATION METHODOLOGY

oil reservoir (Echeverr´ Ciaurri et al., 2011; Bellout et al., 2012) applications. For ıa cases where we can access multiple processors in parallel, we apply Particle Swarm Optimization (PSO), which is a stochastic global search method (Kennedy & Eberhart, 1995). In this chapter, we will also describe the Generalized Pattern Search optimization method, which we will use in Chapter 4 for data assimilation.2.1 Deﬁnition of the optimization problem

**The optimization problem in this work can be expressed generally as:**

where J is the objective function we seek to minimize, u ∈ Ω ⊂ Rn are the optimization variables deﬁning well locations, injection rates and water cycling events (if included), and Ω represents the various bound and inequality constraints deﬁning a convex feasible region for u (see, e.g., Boyd & Vandenberghe, 2004, for the deﬁnition of a convex feasible region).

In many of our optimizations, the objective function J is taken to be the fraction of injected CO2 in the mobile phase after 1000 years (CO2 is injected for 30 years and the system equilibrates for an additional 970 years). Optimizing this objective facilitates safe long-term trapping of CO2. In some examples we will also consider diﬀerent deﬁnitions for J, which will be indicated by an appropriate subscript. These include the mobile CO2 fraction as a time-average over the 1000-year period (Jave ), or over a shorter time frame of 100 years (J100 ). A slightly more complicated deﬁnition (used extensively in Chapter 4), which may be a more accurate indicator of the risk of leakage in some cases, is the time-average of the total mobility of CO2 that reaches the cap rock (Jλ ). Mathematically, this objective is expressed as

where ρg is the gas (CO2 ) phase density, krg is the gas relative permeability (which increases more than linearly with gas saturation Sg ), and µg is the gas viscosity.

Other choices for the objective function are also possible. For example, we might seek to minimize the fraction of CO2 in the top layer of the model, or maximize dissolved CO2 after the equilibration period. Other eﬀects, such as average pressure or pressure at the top of the aquifer, could also be included in the objective function.

In real cases, it may also be desirable for the CO2 plume to avoid certain areas of the formation which may have a greater risk of leakage.

The optimization variables u are divided into ﬁve categories, u = [x, r, s, l, p].

Here x describes the locations of N CO2 injection wells, r describes the CO2 injection rates of each well over M time intervals, s denotes the starting times for K brine cycling events, l denotes the duration of each cycling event, and p denotes the proportion of the total injected or produced brine at each of the K cycling events in each of the N cycling wells. A general description for each variable is provided below. Table 2.1 contains summary information for the control variables used in our optimizations.

In order to limit the dimensions of the search space, we restrict the CO2 injectors to be horizontal wells that lie in the bottom layer of the aquifer and extend in either

## 20 CHAPTER 2. OPTIMIZATION METHODOLOGY

the x or y-direction (the orientation of each well is speciﬁed). This means that each injection well is deﬁned by three real variables, denoted (x1, x2, x3 )|w∈{1,...,N } for w w w a total of 3N optimization variables. The ﬁrst two variables, (xw, x2 ), correspond w to the (x, y) location of the heel (starting point) of the well, and x3 speciﬁes the w length of the well. Because wells are prescribed to start and end in the centers of discrete grid blocks, these optimization parameters are treated as integer values. We accomplish this in HJDS by controlling search distances, as discussed in Kolda et al.(2003), and in PSO by rounding to the next integer value using the ceiling function.

The placement variables are constrained to ensure that wells fall within the allowable (core) region of the aquifer and that they are neither too long or too short. More general wells (e.g., deviated or multilateral), as described in Yeten et al. (2003) and Onwunalu & Durlofsky (2010), could be considered, though this would increase the computational requirements of the optimization problem.

The injection rate optimization variables r represent fractions of the total injection rate and are denoted rw,t for w ∈ {1,..., N − 1} and t ∈ {1,..., M }, giving a total of M (N − 1) independent control variables. Note that the rate fractions for well N can be omitted as control variables, because they can be written as linear combinations of the other rate fractions, such that

Furthermore, injection rate fractions are bounded on the interval [0, rmax ]. This bound must also be enforced for well N, which gives rise to the linear inequality constraints

where rmax is the maximum fraction of the total injected CO2 associated with any individual well (in our optimizations we set rmax = 0.4).

In our modeling, brine cycling events entail injecting brine in the top layer of the model while simultaneously producing the same volume of brine from the bottom of the aquifer. Overproduction relative to injection could be applied for pressure management and the pressure level could be formally optimized, though this

## 2.1. DEFINITION OF THE OPTIMIZATION PROBLEM

was not considered here. Our numerical investigations showed that brine cycling is most eﬀective when performed as separate (distinct) events. Thus, the optimization variables for brine cycling deﬁne the timing, duration and the pore volume fraction produced/injected at each well for K cycling events. Note that cycling volumes are expressed as fractions of the pore volume of the ‘core’ aquifer region, deﬁned in Section 2.3.Again, we introduce some assumptions to limit the dimensions of the search space.

Namely, we specify that there are N brine injectors, which are in the top layer of the model and are aligned (in terms of x − y location) with the N CO2 injectors. For each brine injector there is a corresponding brine producer. This conﬁguration shares some similarities with the engineered injection strategy proposed by Anchliya et al.

(2012), shown in Figure 1.3. In our model, however, the brine producers are located in the bottom layer just outside the four corners of the core aquifer region. Thus they do not produce CO2. Figure 2.1 shows the general locations of all wells for cases with brine cycling. We reiterate that the locations of the CO2 injectors are optimization variables, and these locations deﬁne the placement of the brine injectors. Thus we can view the optimization as providing the locations of these well pairs. The production rates and injection rates are speciﬁed to be equal for each brine injector/producer pair.

Given these assumptions, we can characterize the timing, duration and pumping rates for K brine cycling events using K(2 + N ) − 1 optimization variables.

** Figure 2.1: Illustration of well positions for cases with brine cycling **

**22 CHAPTER 2.**

OPTIMIZATION METHODOLOGY

Speciﬁcally, the optimization variables for brine cycling are speciﬁed as: s1,..., sK, denoting the start time of each event, l1,..., lK denoting the duration of each event, and pw,k for w ∈ {1,... N } and k ∈ {1,..., K}, denoting the injection/production pore volume fraction for each well pair during each cycling event. The optimization variable pw,k is expressed as a fraction of the total (speciﬁed) injected pore volume over all cycling events. Consequently, we can remove pN,K as an independent optimization variable by applying the constraintFurthermore, the actual injection/production rates associated with the fractions pw,k are bounded by a maximum value, designated qmax. Here we used qmax = 5 × 104 bbl/day = 7950 m3 /day (this constraint was reached only once in our optimizations).

**This introduces N K convex constraints as follows:**

pw,k 0≤ ≤ qmax, ∀ w ∈ {1,..., N }, k ∈ {1,..., K}. (2.8) lk Note that for w = N and k = K in Equation 2.8, we substitute for pN,K using Equation 2.7.

A number of optimization algorithms could be applied to solve the problem described in Equation 2.1. Because we are using a commercial simulator that does not provide gradients from an adjoint procedure, gradient-based techniques cannot be applied eﬃciently. We thus apply noninvasive gradient-free optimization approaches.

Several such techniques were considered in Echeverr´ Ciaurri et al. (2011) for well ıa control optimization in the context of oil reservoir management. There it was shown that Hooke-Jeeves Direct Search (HJDS) is an appropriate method to use in cases when parallel computational resources are limited. We note that, in practice, parallel resources may be limited due to core availability or by limitations on the number of

2.2. OPTIMIZATION METHODS commercial software licenses. If multiple cores and licenses are available, other (readily parallelizable) direct search methods such as Generalized Pattern Search (GPS) or Mesh Adaptive Direct Search (MADS), or global stochastic search methods such as Genetic Algorithms (GA) or Particle Swarm Optimization (PSO), could also be applied.

2.2 Optimization methods 2.2.1 Hooke-Jeeves Direct Search In this work, our optimizations involve a total of n = 3N + M (N − 1) optimization variables (without brine cycling), or n = 3N + M (N − 1) + 2K + N K − 1 optimization variables (with brine cycling). We consider examples with N = 4, M = 4 or 6, and K = 3 (in cases with brine cycling). Thus, n ranges from 24 to 47 variables.

Problems of this size can be readily addressed using direct search techniques, where the search involves the evaluation of a cost function using a stencil-based approach.

HJDS entails two types of moves through the search space, which are referred to as exploration moves and pattern moves. Search directions, designated di, are aligned with the search coordinates for exploratory moves, though for pattern moves they are in general not aligned with search coordinates. The magnitudes of the components of di deﬁne the step size. Constraints are handled by projecting infeasible points onto the feasible region using the approach described by Kolda et al. (2003).

An illustration of the algorithm, for a minimization in R2, is shown in Figure 2.2.

HJDS ﬁrst performs exploration (stencil-based) moves from an initial point, u0. The search directions are polled sequentially and opportunistically. Speciﬁcally, the ﬁrst exploration computation involves a function evaluation at u0 + d1. If J(u0 + d1 ) J(u0 ), then the move is successful and the next point evaluated is u0 + d1 + d2. If J(u0 + d1 ) J(u0 ), then u0 − d1 is considered. The opportunistic aspect of the search means that, as soon as improvement is achieved in any search direction, we proceed with the next search coordinate. If no improvement in the cost function is found after the exploratory step, then the step size is reduced. Otherwise, once all di, i = 1,..., n are considered, HJDS performs a pattern move. The direction for the pattern move is based on the outcomes of the exploratory moves. Speciﬁcally,

## 24 CHAPTER 2. OPTIMIZATION METHODOLOGY

Figure 2.2: Illustration of exploratory and pattern moves in HJDS. The star represents the optimum (from Echeverr´ Ciaurri et al., 2011) ıa if we designate the location in search space after the exploratory moves as u1, then the next point evaluated in the pattern move is u2 = u0 + 2(u1 − u0 ), as illustrated in Figure 2.2. If J(u2 ) J(u1 ), exploratory moves centered around u2 are then performed. If J(u2 ) J(u1 ), then u2 is discarded and a new polling is performed around u1. If there is still no decrease in J, then the step size is reduced, and polling around u1 is again performed.The algorithm requires an initial guess u0, which we take to be the heuristic base case. Under reasonable assumptions, HJDS can be proven to converge to a local minimum (Kolda et al., 2003). During the ﬁrst stages of the optimization, we use a relatively large step size, which is expected to increase the robustness of the algorithm with respect to noisy objective functions and to enable some amount of global exploration. For the well placement variables, which are required to be integer values for ﬂow simulation (ensuring that they correspond to grid block indices), the minimum perturbation size or tolerance is one grid block. When a variable reaches its minimum tolerance, the step size for that variable remains unchanged for the rest of the optimization. When all search variables have reached their tolerances, and the objective function can no longer be improved, the algorithm is considered to have converged.

Because HJDS only achieves a local minimum, we generally run our optimizations

2.2. OPTIMIZATION METHODS several times using diﬀerent u0. It is also worth noting that HJDS can be readily hybridized with stochastic search algorithms, such as PSO, to incorporate a global search component. See, for example, Aliyev (2011) or Isebor et al. (2013) for a PSO-MADS hybrid. Hybrid approaches are not applied here, though this would be straightforward to investigate. We found HJDS to be an eﬀective technique for optimization when parallel resources could not be used (in our work, this was the case due to limited licenses for commercial software). When multiple processors and licenses were available, however, we found PSO to be similarly eﬀective and to run in signiﬁcantly less elapsed time. We now describe the PSO algorithm used in this work.