Bayesian Optimization for Field-Scale Geological Carbon Storage

Xueying Lu*,Kirk E.Jordn,Mry F.WheelerEdwrd O.Pyzer-Knpp,Mtthew Bentn

a Center for Subsurface Modeling, Oden Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX 78712, USA

b IBM Research United Kingdom, Warrington WA4 4AD, UK

Keywords:Compositional flow Bayesian optimization Geological carbon storage CCUS Machine learning AI for science

A B S T R A C T We present a framework that couples a high-fidelity compositional reservoir simulator with Bayesian optimization (BO) for injection well scheduling optimization in geological carbon sequestration. This work represents one of the first at tempts to apply BO and high-fidelity physics models to geological carbon storage. The implicit parallel accurate reservoir simulator (IPARS) is utilized to accurately capture the underlying physical processes during CO2 sequestration. IPARS provides a framework for several flow and mechanics models and thus supports both stand-alone and coupled simulations.In this work, we use the compositional flow module to simulate the geological carbon storage process.The compositional flow model, which includes a hysteretic three-phase relative permeability model,accounts for three major CO2 trapping mechanisms: structural trapping, residual gas trapping, and solubility trapping. Furthermore, IPARS is coupled to the International Business Machines (IBM)Corporation Bayesian Optimization Accelerator (BOA) for parallel optimizations of CO2 injection strategies during field-scale CO2 sequestration. BO builds a probabilistic surrogate for the objective function using a Bayesian machine learning algorithm—the Gaussian process regression, and then uses an acquisition function that leverages the uncertainty in the surrogate to decide where to sample. The IBM BOA addresses the three weaknesses of standard BO that limits its scalability in that IBM BOA supports parallel (batch) executions, scales better for high-dimensional problems, and is more robust to initializations. We demonstrate these merits by applying the algorithm in the optimization of the CO2 injection schedule in the Cranfield site in Mississippi, USA, using field data. The optimized injection schedule achieves 16% more gas storage volume and 56% less water/surfactant usage compared with the baseline. The performance of BO is compared with that of a genetic algorithm (GA) and a covariance matrix adaptation (CMA)-evolution strategy (ES). The results demonstrate the superior performance of BO, in that it achieves a competitive objective function value with over 60% fewer forward model evaluations.

1. Introduction

Under the framework of the Paris Agreement,achieving carbon neutrality by the middle of the century is the fundamental solution to cope with the climate crisis. Along with electrification, hydrogen, and sustainable bioenergy, carbon capture, utilization, and storage(CCUS)is key in achieving a net-zero energy system.CCUS is the only group of technology that contributes both to removing hard-to-abate CO2emissions and directly reducing emissions in the key sectors [1]. In geological sequestration, CO2is captured from a power plant or industrial facility and compressed;it is then transported to and injected as a supercritical fluid into depleted oil and gas reservoirs and deep saline aquifers for long-term storage.The basis for CCUS’s potential is the huge global storage capacity that exists in geological formations and the availability and proximity of potential injection sites to power generation plants. However, such injections present significant technical challenges in regard to ensuring the safety of injection operations and minimizing the probability of leakage through geological faults and naturally connected cracks on a time scale of hundreds or even thousands of years, with controllable costs. In addressing these challenges, accurate prediction of the behavior of injected CO2is important to the long-term success of carbon sequestration,because even small leakage rates over long time periods can unravel the positive outcomes of net sequestered CO2[2-5].Furthermore,reservoirs’ geomechanical responses associated with the impact of high rates of injection of CO2on natural fractures are crucial in determining earthquake effects[6,7].Here,accurate prediction of the fate of injected CO2and the related responses of geological formations must encompass physical models governed by multiphase flow,multicomponent transport,rock mechanics,thermodynamic phase behavior, chemical reactions within both fluid and rock, and the coupling of all these phenomena over multiple temporal and spatial scales. This effort requires high accuracy in the physical models used and in their corresponding numerical approximations. Aside from geological carbon storage, CO2enhanced oil recovery (EOR), in which oil recovery is improved by injecting CO2into oil reservoirs,is another key method in CCUS.Recent work[8]has demonstrated that the amount of sequestrated CO2in storage-driven CO2EOR can exceed the amount of emissions from burning the produced oil,resulting in net-zero or even negative CO2emissions.

In this work, we use the implicit parallel accurate reservoir simulator (IPARS) developed at the Center for Subsurface Modeling at the University of Texas at Austin to simulate the geological carbon sequestration processes. IPARS provides a framework that supports both stand-alone and coupled simulations of several physical models. The models include single and two-phase flow, a black oil model, an equation-of-state (EOS)compositional flow model, a thermal and geomechanics model,and so forth [9-13]. IPARS has been used extensively for field-scale benchmark and experimental studies related to CO2storage in geological formulations [13-16]. More specifically,the EOS compositional flow model is adopted to simulate the CO2storage process in a deep saline aquifer. Both CO2and brine are modeled as hydrocarbon components in the compositional model to account for the effect of CO2dissolution into brine. The Peng-Robinson EOS [17] is used for CO2-brine phase behavior and property calculations. A flash calculation follows to determine the mole fractions of CO2and water in two equilibrium phases. In particular, this module is enhanced to model CO2-brine phase behavior for CO2geological sequestration applications [18,19]. The compositional flow module contains advanced petrophysical models, including a novel three-phase,hysteretic relative permeability model and an implicit texture foam model for the accurate characterization of cyclic injection processes and foam-assisted recovery processes [5,20,21]. These models account for the residual gas-trapping effect during geological carbon storage processes. IPARS has also been used for coupled flow and mechanics simulations to study poro-elastic and poro-plastic effects during geological carbon storage[5,22]. The flow modules in IPARS are formulated on general hexahedral grids using the mass-conservative multipoint flux mixed finite-element (MFMFE) method [10]. Moreover, IPARS supports multiple preconditioners and nonlinear solvers, such as Bi-conjugate gradient stabilized (BCGS) method, algebraic multigrid (AMG), and generalized minimal residual (GMRES)method. The simulator can handle millions of grid blocks with a message passing interface (MPI) library for parallel distributed memory computations in a high-performance computing environment [11].

Design and control problems in reservoir well management are challenging, as the forward problem involves coupled multiphysics simulations, which are highly nonlinear, multiscale, computationally expensive, and under geological uncertainties. From an optimization point of view, derivative-free optimizations are popular. These methods include particle swarm, simulated annealing (SA), genetic algorithms (GAs), and evolution strategy(ES), to name a few [23-25]. These global optimization methods enjoy the flexibility of using existing reservoir simulation software; hence, they have been extensively utilized for case studies. Another group of methods utilize stochastic approximations of the gradients, such as simultaneous perturbation stochastic approximation (SPSA) [26] and the stochastic simplex approximate gradient (StoSAG) [27], to name a few. Adjoint methods have also been applied to reservoir well optimizations[28,29]. Nevertheless, large-scale applications are still computationally demanding.

Recently,machine learning,reduced-order modeling,and other types of surrogate models have also been extensively integrated into subsurface simulations for speedups [30-32]. Bayesian optimization (BO) has emerged as a powerful solution for control and optimization problems when the objective functions are expensive to evaluate or potentially intractable. BO has been shown to be successful in machine learning for hyperparameter tuning [33], as well as in several scientific domains such as material design,robotics,and environmental monitoring[34-36],especially for low-dimensional problems [37]. Compared with other derivative-free methods, BO leverages a machine learning technique, Gaussian process regression, which permits analytical tractability when performing Bayesian inference.The combination of an acquisition function to statistical inference allows an exploration-exploitation trade-off in searching the surrogate solution space. BO has been used for uncertainty quantification for a benchmark reservoir model [38]. In this work, we adopt the BO framework to optimize the injection well scheduling in geological carbon sequestration with data from the Cranfield reservoir in Mississippi, USA.

This paper is organized as follows: Section 2 introduces the mathematical modeling of carbon sequestration processes, Section 3 introduces BO and its coupling with the reservoir simulation,the experiments and results are discussed in Section 4, and Section 5 presents the conclusions.

2. Mathematical modeling of the carbon sequestration process

2.1. Compositional flow model

In this work,we adopt the parallel,EOS,and compositional flow module in IPARS to model the process of geological carbon storage in a depleted reservoir or a deep saline aquifer. We assume that isothermal flow, water, and Nchydrocarbon components form a three-phase flow system—namely, an aqueous phase, a nonaqueous liquid phase, and a gaseous phase. Let φ denote the current fluid fraction (i.e., porosity), Sα the saturation of phase α, ραthe mass density of phase α, qiα the injection/production rate of component i in phase α, ξiα the mole fraction of component i in phase α, uα the volumetric velocity of phase α, and Diα the diffusion-dispersion tensor of component i in phase α. The mass conservation equation of component i for a multiphase flow system is then given as follows:

where t is the time variable;o is oil phase;w is water phase;g is gas phase.

Let κ denote the absolute permeability tensor of the porous rock matrix, krαthe relative permeability,μαthe viscosity, and pαthe reference phase pressure; g is the magnitude of the gravitational acceleration and z is the depth. The volumetric velocity of phase α, uα, is described by Darcy’s law:

The capillary pressure equations,auxiliary equations,and initial and boundary conditions are supplemented to close the system[39,40]. The primary variables chosen to solve the system are a reference phase pressure pαand the Nccomponent concentrations S1， S2， ...， SNc. A slightly compressible and a cubic Peng-Robinson EOS [17] are used for the water and hydrocarbon phases, respectively.The phase equilibrium is calculated using the Rachford-Rice equation and the isofugacity criteria.

2.2. Hysteretic relative permeability and foam models

One of the lessons learned from gas injection in oil reservoirs is that the gas volumetric swept efficiency is inherently low due to gravity override caused by low gas viscosity, density, and formation heterogeneity. Therefore, methods to improve the volumetric sweep efficiency are valuable. Gas mobility can be reduced by the alternative injection of gas and water (WAG). During the injection of a water slug, gas is trapped by capillary forces in high gas saturation regions,so that gas flow is diverted into regions with lower gas saturation in the subsequent gas slug [41]. However, WAG effectiveness is limited, and WAG is not very efficient in highly heterogeneous or thick reservoirs. Hence, surfactant-alternatinggas (SAG) can be introduced as a second remedy to reduce gas mobility.The in situ generated foam increases the gas apparent viscosity and blocks high permeability streaks, thereby diverting the flow from high permeability to low permeability and unswept regions [42].

The UTHYST model [20] offers a simple approach to calculate the cycle-dependent irreversible hysteresis behavior of the relative permeability for any phase.It models the monotonically increasing trapped saturation including caplillary-desaturation effect at high trapping numbers, also uses a dynamic Land coefficient in threephase porous media flow. More details on the model and validations can be found in Ref. [20].

An improved foam model based on the implicit texture foam model in Computer Modelling Group (CMG) STARS has recently been developed and implemented in IPARS [21]. The model proposed a new foam generation function that captures ①the transition of coarse foam to strong foam, ②the foam shear-thinning rheology, and ③ the foam generation hysteresis observed through experiments. The model has been calibrated and validated with experimental data. More details can be found in Ref. [21].

3. Bayesian optimization

Consider a general optimization problem:

where Γ（·） and Hν（·） are the Gamma function and the Bessel function of order ν, respectively.

Not a house was to be seen; the only shelter he could get was the hollow trunk of a great tree, and there he crouched19 all the night which seemed to him the longest he had ever known

With the Gaussian process representing the belief about the objective function f, an acquisition function, denoted by a:X →R+,is used to select what point in X will be evaluated next via a proxy optimization:

where xbestis the best candidate among x1， x2， ...， xn; σn（x ） is the posterior standard deviation of the Gaussian process at n iterations;φ（·） is the probability density function of the standard normal distribution; Φ（·） is the cumulative distribution function of the standard normal distribution. EI encourages both exploration and exploitation—intuitively, it takes higher values when the surrogate mean is high (exploitation) and the surrogate variance is high (exploration). The hyperparameter ε ＞0 allows more flexibility in the exploration-exploitation trade-off [46].

It is shown in Ref. [47] that CI is much more robust to poor priors.

In summary,the pseudo-code of the BO is given in Algorithm 1[48].

Algorithm 1. BO template [48].1: Initialize the probabilistic surrogate with n0 points,D0 = x1，y1■（ ）， x2，y2■ ■（ ）， ...， xn0，yn0■;2: for n=0 to N do 3: Select new query point xn+1 by optimizing acquisition function a,xn+1 =argmaxxa x;Dn（ ）;4: Evaluate yn+1 =f （xn+1）;5: Augment data Dn+1 =Dn ∪ xn+1，yn+1■■ ■■;6: Update the probabilistic mode.

3.1. Parallel and distributed BO

In this work, we adopt a batch implementation of BO based on parallel Thompson sampling (TS) [49]. The acquisition function of TS can be formulated as follows [50]:

where f is parametrized by θ and f（n）～GP （θ|Dn）.

Algorithm 2. Parallel and distributed BO with TS [49].1: Initialize the probabilistic surrogate with n0 points,D0 = x1，y1■. Batch size B.2: for n=0 to N do 3: Compute current posterior p（θ|Dn）.4: for b=1 to B do 5: Sample a random realization of the objective function from its posterior by sampling θ from p（θ|Dn）6: Select new query point xn+1，b ←argmaxx?Dn fθ （x ）7: Evaluate yn+1，b =f （xn+1，b）8: Augment data Dn+1 =Dn ∪ xn+1，yn+1■（ ）， x2，y2■ ■（ ）， ...， xn0，yn0■■ ■■B b=1

3.2. Coupling with the reservoir simulator

The International Business Machines (IBM) Corporation Bayesian Optimization Accelerator(BOA)is hosted on the IBM server.It is coupled to the reservoir simulator IPARS,which runs on an in-house cluster, via an interface function. Fig. 1 illustrates the schematic of the BOA-IPARS coupling.For each batch,the reservoir simulations are run on a local computer cluster in parallel,concurrently; when the simulations finish, results are post-processed to calculate the objective function values. The objective function values are then passed to the IBM cloud for the BOA computations;next, new evaluation points suggested by BOA are returned to the cluster, and simulations are initiated for the next batch.

4. Numerical experiments

4.1. The Cranfield reservoir model

The Cranfield site, which is located in Mississippi, USA, is a depleted reservoir with a salt-cored dome geological structure.The Bureau of Economic Geology (BEG) research team at the University of Texas initiated a field-scale CO2sequestration pilot project at the depleted Cranfield site in December 2009 [15]. By 2015,around 0.5 million metric tons of CO2had been injected into the Cranfield field. The reservoir model contains an 80 ft × 7200 ft × 7200 ft (1 ft = 0.3048 m) domain at a depth of 9901 ft, which is discretized into a distorted hexahedral mesh of 20 × 144 × 144 elements with block dimensions of 4 ft × 50 ft × 50 ft. Five injection wells are simulated. Since the flow model entails no-flow boundary conditions, eight pseudo-production wells are assigned to the boundaries of the reservoir to mimic open-boundary conditions (Fig. 2). The pseudo-production wells are constrained at a constant bottom-hole pressure equal to the initial reservoir pressure, which is 4650 psi (1 psi = 6.894757 kPa); the initial reservoir temperature is 125 °C. The reservoir is assumed to be initially saturated with brine. Both CO2and brine are modeled as hydrocarbon components in the compositional model to account for the effect of CO2dissolution into the brine.The Peng-Robinson EOS [17] is used for the CO2-brine phase behavior and the property calculations. A flash calculation follows to determine the mole fractions of CO2and water in the two equilibrium phases.

Table 1 [15] provides a summary of the pressure, volume, and temperature (PVT) data. The binary interaction coefficients and volume shift parameters used in the EOS are modified according to published correlations that were fitted on experimental data for the solubility of CO2in brine and the density of brine[18,19,51].The history-matched full-tensor permeability,porosity,relative permeability, and capillary pressure curves are used [15](Figs. 3-5). The Cranfield reservoir data has been used for many numerical simulation studies on geological carbon storage[5,13,15,22].

4.2. Baseline simulations

Fig. 1. Schematic of the BOA-IPARS coupling. nbatch: the number of simulation jobs in a batch.

Fig. 2. Distribution of the wells in the reservoir model, the five injectors are modeled from real injection wells in the field,the pseudo boundary wells mimic the open flow boundary conditions on three boundaries of the reservoir model.

The effect of gas mobility control techniques, including WAG and SAG [25], is demonstrated by simulating the four injection schemes summarized in Table 2: continuous CO2injection, WAG without relative permeability hysteresis modeling, WAG with relative permeability hysteresis modeling, and SAG with relative permeability hysteresis modeling. In the four simulation studies,the gas injections are rate-controlled, such that the total amount of CO2injected in the four cases is equal. The design allows an unbiased comparison of CO2storage volume for each injection and simulation strategy. The simulation results for the CO2distribution at the end of the injections—both views from the top and the bottom—are presented in Fig. 6 for each of the four cases. As expected, during the continuous CO2injection (Fig. 6(a)), the gas plume migrates upwards due to low gas density and viscosity until it reaches the sealing caprock;it then migrates horizontally,resulting in poor sweep efficiency at the bottom part of the reservoir.This scenario is unfavorable from a sequestration standpoint:After 16 years of injection, most of the injected gas is produced through the boundary wells. The significant effect of relative permeability hysteresis modeling is illustrated by comparing the performances of case (Fig. 6(b)) and case (Fig. 6(c)). Both cases simulate a 20-year WAG process; case (Fig. 6(c)) enables relative permeability hysteresis modeling,while case(Fig.6(b))does not.The prediction in case(Fig.6(c))shows a broader CO2swept area at the bottom of the reservoir in comparison with case (Fig. 6(b)). The relative permeability hysteresis accounts for the capillary trapping of the gas phase, leading to the prediction of more CO2storage volume. As summarized in Table 3,without including the relative permeability effect, the WAG process results in 9% less cumulative CO2storage volume compared with continuous CO2injection; when the relative permeability effect is considered, the WAG process can result in a 55% increase in the cumulative CO2storage volume.

The prominent effect of SAG in gas mobility control is shown in Fig. 6(d). As the liquid/gas mobility ratio is reduced by the in situ generation of foam, the gas phase can invade the lowpermeability area and overcome gravity segregation. The vertical gas sweep efficiency is significantly improved in comparison with both the continuous CO2injection and WAG. The amount of CO2storage volume in each of the injection cases is summarized in Table 3, showing that SAG significantly improves the CO2storage volume by 77% compared with continuous CO2injection. The SAG simulation with the injection schedule in Table 2 is set as the baseline for the injection cycle optimization in the next subsection.

Fig. 3. Reservoir initial porosity.

Fig. 4. Reservoir vertical permeability in logarithmic scale (logpermx). md:millidarcy.

Fig. 5. Relative permeability and capillary pressure curves. krg: gas phase relative permeability; krw: water phase relative permeability; Pc: capillary pressure.

Table 1 EOS data for CO2 sequestration in Cranfield [15].

Table 2 Summary of CO2 injection schemes.

Fig. 6. Simulation results of CO2 distribution as gas saturation (Sgas) at the end of injection for each of the injection cases. (a) Continuous CO2 injection; (b) WAG without relative permeability hysteresis; (c) WAG with relative permeability hysteresis; (d) SAG with relative permeability hysteresis.

Table 3 Summary of CO2 storage volumes for different injection schemes.

4.2. BO for gas/surfactant injection scheduling

The objective is to optimize the gas/surfactant injection cycles in order to maximize the cumulative carbon storage volume in the reservoir model described above within a specific injection period.Let T denote the total injection time,rinj，ithe CO2injection rate of the injection well i, Ninjthe total number of the injection wells,rprod，ithe CO2production rate of the production well i,and Nprodthe total number of production wells. The optimization problem is then formulated as follows:

In the BO experiments, we used the Matérn 52 covariance kernel for the Gaussian process and the CI acquisition function. A number of 15 function evaluations were assigned to initialize the Gaussian process.Aside from the initialization,the computing budget was set to a maximum number of 200 function evaluations.For the batch BO, we used a batch size of five. We also experimented with GAs and the covariance matrix adaptation (CMA)-evolution strategy (ES) [25] on the same optimization problem; for the implementation of these algorithms, we used the open-source ES package Distributed Evolutionary Algorithms in Python (DEAP)[52].

The results from the best experiment of BO are summarized in Table 4.It is noticeable that the baseline SAG scheme is targeted at five injection cycles, starting with a gas injection cycle. The optimized SAG scheme only has two injection cycles, starting with a surfactant injection cycle.The results are consistent with the physical mechanism of in situ foam generation—an initial cycle with surfactant injection will facilitate foam generation at an early stage of the entire injection scheme. Compared with the baseline SAG schedule, the optimized schedule achieves 16% more CO2storage volume within a 20-year injection plan and reduces 56% of the water and surfactant usage.

We test the robustness of the BO algorithm by comparing its performance with the other two algorithms in repeated experiments. The experiments with each algorithm are repeated three times. Fig. 7 compares the performance of the best runs among the repeated experiments for each of the three algorithms. The experiments show that BO finds a better solution than the GA and CMA-ES,with a smaller number of function evaluations.More specifically, BO reaches its best observation, 1.31 × 108Mscf (1 Mscf = 28.317 m3), after 73 function evaluations, whereas CMAES reaches 1.28 × 108Mscf after 183 function evaluations. BO achieves competitive objective function values with 60% fewer function evaluations. These results demonstrate the efficiency of BO. Fig. 8 plots the overall performance of the three algorithms through repeated experiments. The solid line represents the mean of the function values, and the shaded area represents one standard deviation. These repeated experiments demonstrate the robustness of the algorithmic advantages of the batch BO.

Fig. 7. Comparison of the performance of the best runs by BO, GA, and CMA-ES.

Fig. 8. Comparison of the overall performance of BO, GA, and CMA-ES.

Table 4 Optimized gas/surfactant schemes via BO compared with continuous gas injection and the baseline injection scheme.

5. Conclusions and discussion

The proposed framework represents a first attempt to incorporate high-fidelity physical models and machine learning techniques that scale up to field-scale applications of geological carbon sequestration. We coupled a high-fidelity compositional reservoir simulator,IPARS,with the IBM BOA to optimize injection well scheduling in order to maximize the cumulative CO2storage volume with in situ generated foam as a gas mobility control technique. In this framework, the forward simulations using IPARS strictly honor the complicated physics during geological carbon sequestration. The compositional flow model, which includes a hysteretic, relative permeability model and a foam model, can simulate the three major trapping mechanisms of CO2storage—namely, structural trapping, residual trapping, and solubility trapping. BO was demonstrated to be successful in the design and to control where the objective functions are expensive to evaluate.Here, the BO consists of two components: a Gaussian process machine learning model that builds a Bayesian surrogate for the objective function, and an acquisition function that decides the sequence of sampling points. The Gaussian process allows the usage of Bayesian statistics for analysis and grants tractability regarding the uncertainties in the unknown objective function values where evaluations are not available, which is used by the acquisition function to balance the exploration-exploitation trade-off. The combination of these two components enables an efficient search of the surrogate solution space.

We applied this framework to field-scale carbon storage with field data from the Cranfield reservoir site to design injection cycles for foam-assisted carbon storage. Our numerical experiments showed that the optimized injection scheme achieves 16%more gas storage volume and 56% less water/surfactant usage compared with the baseline case.We further conducted a comparison study to examine the performance and robustness of BO and compare them with those of other commonly used algorithms,including a GA and CMA-ES. The results demonstrate the superior efficiency and robustness of the BO algorithm, in that it achieves competitive objective function values with a much smaller number of function evaluations. The results show the promise of further applications of BO for other control and well management applications in CCUS projects, such as optimization with economic constraints, optimization for well placements, and so forth. Since IPARS includes geomechanics modules that are efficiently coupled with the compositional flow model, the framework can be easily extended to include geomechanics effects during carbon storage.Furthermore,the framework can be extended to consider BO under geological uncertainties.

Acknowledgments

Xueying Lu and Mary F. Wheeler are have been supported under the Center for Subsurface Modeling Affiliates Program,United States of America and the National Science Foundation,United States of America (1911320, Collaborative Research: High-Fidelity Modeling of Poromechanics with Strong Discontinuities).