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

The optimization determines the locations of four horizontal CO2 injection wells, each of maximum length 1304 m (the length of 3 grid blocks). The maximum amount of CO2 that a well can inject in any time interval is 40% of the total CO2 injected (this maximum acts to limit the pressure increase in the vicinity of the well). Four horizontal water injectors are also used along with producers in boundary blocks in cases with water cycling events, as described previously (see Figure 2.1).

** Figure 2.5: 39 × 39 × 8 areal grid with 25 × 25 × 8 storage aquifer shown as white square in the center The aquifer model is heterogeneous with porosity φ generated using sequential Gaussian simulation (Remy et al.**

, 2009). Porosity is normally distributed with mean ¯ φ = 0.2 and standard deviation σφ = 0.025. The porosity distribution is truncated on the interval [0.15, 0.25]. A spherical variogram model is used, with correlation lengths of cx = 8∆x (in Chapter 3) and cx = 24∆x (in Chapters 4 and 5) in the x-direction, cy = 8∆y in the y-direction, and cz = ∆z in the z-direction, where ∆x,

2.3. AQUIFER MODEL Figure 2.6: 25 × 25 × 8 storage aquifer showing log-permeability ﬁeld. Injection wells can be placed within the inscribed 15 × 15 × 8 ‘core’ region ∆y, ∆z designate grid-block dimensions. This deﬁnes the prior realizations in cases with geological uncertainty. When injection wells and sensor wells are drilled, the porosity values from a particular realization, taken to be the actual (true) model, are included as hard data.

Permeability (kx ) is related to porosity using the relation ¯ φ−φ kx = exp a + b, (2.12) σφ where a = 4.2 and b = 1.5 in the model in Chapter 3, resulting in a log-normal permeability distribution for kx with a median kx of 67 md and minimum and maximum values of 3 md and 1300 md. Examples in Chapters 4 and 5 use a slightly more heterogeneous model, with a = 4.5 and b = 2 in Equation 2.12. This results in a log-normal permeability distribution between about 2 md and 5000 md with a median of 90 md. In all examples we specify ky = kx and kz = 0.1kx.

Relative permeability curves, shown in Figure 2.7, were derived by applying the Brooks-Corey parameters found by Krevor et al. (2012) for a Berea sandstone. The connate water saturation is Swc = 0.11 with corresponding maximum CO2 relative permeability of krg (Swc ) = 0.95. The Brooks-Corey exponents are Nw = 6 and Ng = 5. Note that krw is very low for Sw ≤ 0.5 (i.e., Sg ≥ 0.5). Hysteresis is represented using Killough’s method (Killough, 1976) with a Lands trapping coeﬃcient

## 34 CHAPTER 2. OPTIMIZATION METHODOLOGY

** Figure 2.7: Relative permeability curves (Krevor et al.**

, 2012) for water (drainage only) and CO2 (drainage and imbibition) of 1. Capillary pressure eﬀects are not included in our models. Limited investigation showed that the impact of capillary pressure was small if the capillary pressure curves were taken to be the same in all grid blocks. However, it has been shown in Saadatpoor et al. (2010) and Krause et al. (2009) that heterogeneity in capillary pressure can be an important eﬀect in CO2 trapping, and this should be considered in future work (though we note that capillary heterogeneity can signiﬁcantly increase simulation time, and may also require much ﬁner grid resolution).

All simulations are performed using the CO2STORE option in the ECLIPSE simulator. The equation of state calculations within the simulator match experimental data from Spycher & Pruess (2005). When CO2 dissolves in the aqueous phase, the increased density leads to gravity-driven transport, though this eﬀect may be underestimated due to the large grid blocks used in our models. The initial salinity of the in-situ water is 10,000 ppm. Salt precipitation is not modeled, though water mobility is aﬀected by salt concentration. As we mentioned previously, we do not include mineralization in this study. Our simulations also do not resolve detailed density-driven instabilities, as computed by, e.g., Riaz et al. (2006). Incorporating such eﬀects in

2.3. AQUIFER MODEL practical engineering models could be accomplished in theory by computing appropriate coarse-scale ﬂux functions from detailed ﬁne-scale results, but robust procedures for this upscaling are not yet available.

The simulations in this work are all run fully implicitly, with time steps of about 50 days during the injection period and up to 3650 days, or 10 years, after injection is completed (ﬂuid movement slows dramatically during this time so we can use large time steps). More time steps are required for cases with brine cycling, since this entails additional injection. We typically use approximately 350 time steps for 1000-year simulations without brine cycling and 450 time steps for 1000-year simulations with brine cycling. Run times for a single 1000-year simulation are around 3-5 minutes, and up to 7 minutes for cases with brine cycling. In some optimization examples, we run 500-year or 100-year simulations. These typically reduce the run time of a single simulation by about 30 seconds to one minute. Some data assimilation examples require only 10-year, 30-year or 100-year simulations, and none of these include brine cycling. Run times for those examples are typically around 1 minute, 2 minutes and 3 minutes, respectively.

36 CHAPTER 2. OPTIMIZATION METHODOLOGY Chapter 3

**Optimization results with knowngeology**

In this chapter we present extensive optimization results for optimization examples where the geology is known. The objective function in most cases is the fraction of mobile CO2 in the formation after 1000 years, though other objective functions and shorter time frames are also considered. Results for cases without and with brine cycling are provided. Hooke-Jeeves Direct Search (HJDS) is applied unless otherwise indicated (when this work was performed, the requisite simulations were run serially because we did not have access to multiple software licenses). We also investigate the sensitivity of optimization results to grid reﬁnement.

All examples in this chapter use the geological model depicted in Figure 2.6. This model was produced by Sequential Gaussian Simulation using the Stanford Geostatistical Modeling Software (Remy et al., 2009). As stated previously, porosity values ¯ are normally distributed with mean φ = 0.2, and standard deviation σφ = 0.025.

Permeability is related to porosity via Equation 2.12, with a = 4.2 and b = 1.5. We note that diﬀerent values for a and b will be used for models in Chapters 4 and 5.

3.1 Optimization without brine cycling We ﬁrst consider a case with no water cycling. Because HJDS is a local optimizer, it is useful to run the optimization multiple times. Results for three optimization runs,

## 38 CHAPTER 3. OPTIMIZATION RESULTS WITH KNOWN GEOLOGY

using diﬀerent initial guesses for well placement parameters x0 but the same initial guess for injection rate parameters (rw,t = 0.25, or equal fractions at all time), are shown in Figure 3.1. The ﬁgure shows the best solution achieved up to the current iteration. It is evident that each optimization run gives a diﬀerent result, indicating the presence of multiple local optima. The best initial guess has a mobile CO2 fraction of around 0.30, whereas the best optimized solution, Run 1, has a mobile CO2 fraction of 0.22. Thus we see clear improvement from the optimization.Particle Swarm Optimization (PSO) was also applied for this example and it provided comparable optimal solutions. While PSO required fewer iterations than HJDS (100 compared with around 500), it performed more function evaluations (2000 compared with around 500) because each iteration entailed the simulation of 20 potential solutions (particles). Because PSO, in contrast to HJDS, readily parallelizes, it should be considered if parallel resources (and simulation licenses) are available.

The initial guesses for the well locations for the three runs are depicted in Figure 3.2(a). Each conﬁguration forms a nearly square pattern, with two of the wells extending in the north-south direction and two in the east-west direction. Figure 3.2(b) shows the optimized well locations for each run. Although the solutions are clearly diﬀerent, many of the wells are positioned in the same general areas, suggesting that certain injection locations lead to reduced mobile CO2. This is presumably related to the detailed structure of the permeability ﬁeld, although it is diﬃcult to identify direct relationships between particular (local) permeability structures and optimal well positions due to the complicated well-to-well interactions and time-varying injection rates. We do, however, expect optimal well locations and injection schedules to lead to signiﬁcant horizontal migration of CO2 and to saturation distributions that facilitate residual and dissolution trapping.

We now further investigate the characteristics of the best solution achieved (the optimized result for Run 1). Figures 3.3(a) and 3.3(b) show the well locations along with the vertically-averaged mobile CO2 after 1000 years for the initial guess and optimized solutions. Note that the wells lie in the core aquifer, and thus honor the optimization constraints, and that the injected CO2 has not proceeded to the outer boundary (which deﬁnes the edge of the central portion of the aquifer). The CO2 injection fractions for these two solutions at the four injection periods are presented

## 3.1. OPTIMIZATION WITHOUT BRINE CYCLING

0.3 0.28 0.26 0.24 0.22

in Figures 3.4(a) and 3.4(b). It is evident from Figure 3.4(b) that Wells 2 and 3 inject considerably more CO2 than Wells 1 and 4. This is also reﬂected in Figure 3.3(b).

Furthermore, we see that Well 1 injects only during the ﬁrst two periods, while Well 4 is inactive in the ﬁrst period and injects relatively little during the second period.

** Figure 3.3: Well locations and vertically-averaged fraction of mobile CO2 at 1000 years for Run 1.**

The inscribed box denotes the core aquifer, representing the feasible region for CO2 injection wells

0.2 0.2 0.1 0.1

** Figures 3.5(a) and 3.**

5(b) show the evolution of CO2 trapping for the initial guess

## 3.1. OPTIMIZATION WITHOUT BRINE CYCLING

and optimal solution in Run 1. The signiﬁcant reduction in mobile CO2 fraction (from0.34 to 0.22) is achieved through an increase of 0.07 in immobile (residually trapped) CO2 fraction and an increase of 0.05 in dissolved CO2 fraction. In both cases, there is a negligible fraction of immobile CO2 during the injection phase (ﬁrst 30 years).

For the base case, the fraction of mobile CO2 decreases until around 300 years, while in the optimized case it continues to decrease throughout the entire time frame.

Both dissolution and residual trapping are enhanced in the optimized case. The enhancement in dissolution trapping likely occurs because the CO2 ‘plume’ is less continuous in the optimized case, which exposes it to a larger volume of brine. Increased residual trapping may result from the fact that there is more trailing-edge area in the optimized case, and it is in these regions that this type of trapping occurs. It is also possible that the CO2 saturation at the trailing edges is higher in the optimized solution (this may be accomplished by the varying injection rates evident in Figure 3.4(b)), which would increase the degree of hysteresis and thus residual trapping. The ‘disconnectedness’ of the mobile CO2 (which facilitates trapping) in the optimized scenario relative to that in the base case is evident in Figure 3.6, where we display mobile CO2 at 30 and 1000 years for the two cases. The ‘isosurfaces’ in the plots depict contours with a mobile CO2 saturation of 0.1 (these plots appear fairly similar if we contour CO2 saturations down to 0.01). It is clear in the ﬁgure that the optimization leads to multiple disconnected CO2 regions rather than a single large region.

We note ﬁnally that, on average, CO2 migrates upward more slowly in the optimized case, and this also provides for more secure storage. This is evident in Figure 3.7, which displays the total mass fraction of CO2 (in any form) reaching the top layer for the two cases. At the end of the injection phase, the CO2 mass fraction in the top layer is around 0.12 for the initial guess solution, but only 0.03 for the optimized solution. The CO2 mass fraction for the optimized solution remains below that for the base case over the entire period. This is most likely the result of the optimal solution exploiting the detailed heterogeneity structure to increase horizontal migration and increase block saturations during drainage. Both of these eﬀects would result in slowing the upward migration of the CO2 plume, as we observe in the optimal solution.

## 42 CHAPTER 3. OPTIMIZATION RESULTS WITH KNOWN GEOLOGY

0.4 0.4 0.2 0.2