Dynamic characterization of geologic CO2 storage aquifers from monitoring data with ensemble Kalman filter

International Journal of Greenhouse Gas Control 81 (2019) 199–215

Contents lists available at ScienceDirect

International Journal of Greenhouse Gas Control journal homepage: www.elsevier.com/locate/ijggc

Dynamic characterization of geologic CO2 storage aquifers from monitoring data with ensemble Kalman filter

T

Wei Ma, Behnam Jafarpour , Joe Qin ⁎

University of Southern California, USA

ARTICLE INFO

ABSTRACT

Keywords: Geologic CO2 storage CO2 storage monitoring Monitoring data assimilation Ensemble Kalman filter

Monitoring the evolution of the CO2 plume during geologic storage is essential for conformance, verification, and risk assessment and mitigation. Monitoring data also play a critical role in characterizing the storage formation and improving the reliability of predictive models. We investigate the feasibility of using the ensemble Kalman filter (EnKF) data assimilation framework to estimate the hydraulic properties of storage formations and to predict the migration of CO2 plume from monitoring measurements, including transient pressure and saturation data at scattered wells and time-lapse seismic data (modeled as vertically-averaged saturation differences in time). To properly account for the uncertainty in the knowledge about saline aquifer properties, the initial ensemble of formation properties is generated based on uncertain statistical model (variogram) parameters. While integration of data from scattered wells provides limited improvement in reducing the uncertainty in the initial ensemble, assimilation of time-lapse seismic measurements (represented by vertically-averaged saturation differences in time) with the EnKF leads to more noticeable uncertainty reduction and reasonable estimates of the general connectivity trends in aquifer hydraulic properties. The estimation and sensitivity analysis results suggest important differences in filter performance during and after CO2 injection. This difference is attributed to the change in flow behavior and the dominant forces before and after injection (pressure versus gravitational forces, respectively). Additionally, when prior model realizations miss essential flow-related elements (e.g., fractures) in an aquifer, the filter provides out-of-range updates, which could be interpreted as a systematic problem in the filter design, in this case possible inconsistency in the prior models.

1. Introduction Geologic CO2 storage (GSC) has been considered as a potential technology for mitigating the impact of CO2 on climate change (Metz et al., 2005; IPCC, 1996). Injection of supercritical CO2 into saline aquifers is expected to trigger processes that lead to trapping and storage of CO2, mainly, dissolution in brine, immobilization in rock pore space, mineral reactions and precipitation, and hydrodynamic trapping under an impermeable caprock (Bachu et al., 1994; IPCC, 1996; Van der Meer, 1992; Gunter et al., 1997; Pruess and Garcia, 2002; Pruess et al., 2003; Ennis-King and Paterson, 2003; Riaz et al., 2006; Juanes et al., 2006; Oldenburg and Rinaldi, 2011). Numerical simulation offers a tool to model and predict these processes and determine the fate of CO2 after injection into saline aquifers and depleted hydrocarbon reservoirs (Lindeberg, 1997; Hattenbach et al., 1998; Jessen et al., 2001; Oldenburg and Benson, 2002; Pruess and Garcia, 2002; Pruess et al., 2003; Nghiem et al., 2004; Doughty and Pruess, 2004; Kumar et al.,

2005). Several studies have focused on strategies for maximizing CO2 storage by mainly enhancing residual and solubility trapping (Anchliya et al., 2012; Hassanzadeh et al., 2009; Leonenkoy and Keith, 2008; Shamshiri and Jafarpour, 2012; Cameron and Durlofsky, 2012; Cihan et al., 2015). However, efforts to predict the fate of CO2 and optimize the injection design and strategies rely on knowledge of the spatial distribution of aquifer properties. Monitoring the CO2 plume migration plays a key role in detecting early signs of leakage and in demonstrating successful containment (Jenkins et al., 2015). Monitoring technologies can also provide important information for constraining aquifer hydraulic properties (i.e., the parameters of the aquifer model) and for improving the prediction of the CO2 plume displacement in saline aquifers. Optimization of the CO2 storage efficiency requires continuous integration of monitoring data to update simulation models and enhance their predictive capability. Several model calibration techniques have been developed in the groundwater and petroleum engineering literature (Hill and

⁎ Corresponding author at: Chemical and Electrical Engineering, Viterbi School of Engineering, University of Southern California, Bloom Walk, HED 313, Los Angeles, CA, 90089, USA. E-mail address: [email protected] (B. Jafarpour).

https://doi.org/10.1016/j.ijggc.2018.10.009 Received 23 April 2018; Received in revised form 21 September 2018; Accepted 26 October 2018 Available online 09 January 2019 1750-5836/ © 2018 Elsevier Ltd. All rights reserved.

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Tiedeman, 2007; Oliver et al., 2008) that can be used for GCS projects. However, adaptation of these techniques to account for the specific conditions in GCS is needed for their successful application. Some of the main characteristics of the GCS problem include (i) strong nonlinearities due to complex reactive transport processes, partitioning of components, fluid mobility contrasts, and gravity-driven displacement behavior; (ii) availability of spatially sparse data (from wells) that can only provide limited information about the plume evolution behavior; and (iii) poor knowledge of the flow properties in saline aquifers that can cast significant uncertainty in initial geologic descriptions. These considerations are important for not only properly formulating and parametrizing the data assimilation problem, but also studying data needs and designing effective monitoring data acquisition systems. In a recent study, Espinet et al. (2013) used monitoring data from a limited number of wells to estimate aquifer flow properties and predict the spatiotemporal evolution of reservoir pressure and CO2 plume migration in a storage formation. The aquifer model in that study consisted of a finite number of zones with distinct permeability values. The TOUGH2 code, which is a numerical simulator for multiphase (gas and aqueous), multicomponent (different carbon species within each phase) flow and transport, was used to simulate the CO2 injection process. The authors combined a response surface analysis approach with a global optimization algorithm “Stochastic RBF’’ to estimate a few parameters in their model (i.e., the permeability of each zone) from pressure and gas saturation monitoring data at sparse well locations. The zonation approach (Jacquard, 1965; Jahns, 1966) reduces the number of unknown parameters and ensures that the a-priori rock facies distribution is not violated during calibration. On the other hand, zonation may not be desirable when the prior description of the zones is not reliable (McLaughlin and Townley, 1996). In addition, zonation can lead to elimination of potentially important local flow-relevant features (e.g., small-scale fractures or shale lenses). Therefore, a higher-resolution (e.g., grid-level) representations of rock properties to describe and estimate the heterogeneity in model parameters is more desirable. However, an increase in model parameter resolution is only justified if it can be supported by data resolution. To monitor an expanding CO2 plume after injection requires techniques that can track a dynamically evolving CO2 front relatively frequently. Dynamic data from fixed well locations is not likely to be an effective monitoring strategy due to limited information content and significant drilling costs. Past studies have demonstrated the importance of time-lapse seismic data for monitoring CO2 migration behavior and gaining insight about the storage performance of aquifers (White, 2013; Jenkins et al., 2015; Man et al., 2016; Burnisona et al., 2017; Roach et al., 2017). However, in addition to its cost, the effectiveness of timelapse seismic technology can be affected by difficulties related to acquisition, geophysical equipment, data processing and interpretation (Man et al., 2016). The seismic response to CO2 saturation changes exhibit nonlinear behavior, with a strong response to small saturations (5–10%) and a weak response to saturation increases beyond the initial range (Roach et al., 2017). Another difficulty in monitoring the CO2 saturation with time-lapse seismic data is discriminating between the seismic response to pressure increases versus changes in CO2 saturation. The resolution of time-lapse seismic data, especially in the vertical direction, is also limited. In general, the thickness of the storage formation relative to the seismic wavelength and the depth of the formation control the expected vertical resolution of time-lapse seismic data. In shallow and thick formations, it is possible to obtain vertical resolution in CO2 saturation front, while in thin layers only and average CO2 saturation may be estimated. Data resolution and quality of seismic data tend to improve when cross-well tomography is considered. Given the cost of drilling and seismic data acquisition, storage sites with multiple aquifers at different depths (with the same areal zones) are more economical. These difficulties necessitate the use of stochastic seismic data interpretation methods that can account for non-uniqueness and uncertainties in the estimated CO2 saturation changes.

Uncertainty is also a key component in geologic description of the storage formation. A major difficulty in characterization of saline aquifers is the lack of prior knowledge about these formations as they have not been studied widely, unlike hydrocarbon reservoirs or shallow groundwater aquifers. Given the significant costs associated with drilling, limited well data will be available to constrain geologic descriptions of the storage formations. The limited state of knowledge about candidate geologic sites for CO2 storage can have a significant implication on problem formulation and construction of prior geologic models. The existing sources of uncertainty in the data and geologic models necessitate the use of probabilistic characterization and prediction frameworks to account for alternative CO2 plume migration scenarios. The Bayesian inversion framework offers an elegant approach for model calibration and uncertainty quantification (Tarantola, 2005). In practical applications, however, due to computational limits, rigorous Bayesian approaches such as Markov chain Monte-Carlo become infeasible, and approximate method must be used. Among practical sampling-based techniques, the ensemble Kalman filter (EnKF) (Evensen, 1994) has gained popularity among the subsurface flow modeling community as it offers several advantages, including gradient-free implementation, the use of simple (first and second order) statistical information, and ensemble formulation that enables error propagation in nonlinear systems. These properties make the EnKF suitable for large-scale monitoring data integration problems with complex state-space (black-box) models (see (Aanonsen et al., 2009) for a review). The first application of the EnKF to subsurface model calibration was proposed by Nævdal et al. (Nævdal et al., 2002). Since then, the method has been successfully applied to large-scale subsurface characterization inverse problems in groundwater hydrology (e.g., (Chen and Zhang, 2006; Vrugt and Robinson, 2007)), petroleum engineering (e.g., (Aanonsen et al., 2009; Jafarpour and McLaughlin, 2009; Emerick and Reynolds, 2013; Chen and Oliver, 2017; Luo et al., 2018)), and geothermal reservoirs (Tarrahi and Jafarpour, 2012; Vogt et al., 2012; Tarrahi et al., 2015), to assimilate different types of dynamic data (e.g., multiphase flowrate, pressure, tracer, and geophysical data) into geologic models of groundwater aquifers and hydrocarbon reservoirs. In this paper, we study the feasibility of applying EnKF data assimilation to characterize heterogeneous geologic CO2 storage formations and their hydraulic rock properties with different combination of monitoring data types, including well pressure and saturation measurements and time-lapse seismic. We examine the importance of lowresolution time-lapse seismic data, in comparison to pressure and CO2 saturation (Sg ) data from scattered monitoring wells, for characterizing aquifer hydraulic/storage properties, and for predicting the evolution and storage of the injected CO2. While time-lapse seismic data is often represented in terms of amplitude difference between baseline and monitor surveys, in this study we use changes in CO2 saturation (at two different times) to represent time-lapse seismic data. However, to account for the limitation in the resolution of time-lapse seismic data, we use noise-corrupted vertically-averaged changes in CO2 saturation for data assimilation. The complex buoyancy-driven migration of CO2 and the heterogeneity in the storage formation create important differences in the spatial and temporal distribution of parameter sensitivities, and are consequential for the design, monitoring, and characterization of GCS projects. Furthermore, for geologic CO2 storage, where the initial ensemble of saline aquifer properties is expected to be highly uncertain, we evaluate the performance of the filter in constraining the distribution of the formation hydraulic properties from very limited well data and low-resolution time-lapse seismic surveys. In addition, we investigate the ability of the EnKF to detect inconsistency in prior models, for instance when natural or induced fractures with complex nonGaussian distributions are present in the formation, but not included in model realizations. The paper is organized as follows. Section 2 contains a brief and general presentation of the EnKF methodology, followed by a 200

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

description of the specific model setup used in our numerical experiments. In Section 3, the geologic model used for CO2 storage is described. Sections 4 and 5 present the numerical experiments and a discussion of the observed results, respectively. The summary and conclusions drawn from this work are provided in Section 6.

predicted observations to the vector of state variables, can be expressed as:

xi = {

T i ,

diT }T

(2) th

where x i refers to the i realization of the augmented state vector, i denotes the dynamical state variables such as pressure head and phase saturation for the ith model realization, i stands for the static parameters such as permeability and porosity, and di is model−predicted measurements such as well pressure (the superscript T refers to the transpose of a matrix). The EnKF implementation consists of a sequence of forecast and update steps. Fig. 1 shows one such sequence for the CO2 storage model that is used in the next section. In the forecast step, each realization of the model parameters is used to predict the corresponding realizations of the state variables at the next time step, following the non-linear state propagation function:

2. Data assimilation with ensemble Kalman filter The core principle behind EnKF is using physical correlations (computed via reservoir simulation) between input model parameters and their corresponding predicted responses to update uncertain model parameters based on observed data from the field. The physical correlations are approximated by (i) generating several initial models with different (uncertain) input parameters, (ii) using each parameter realization with reservoir simulation to estimate the corresponding (ensemble of) response data, and (iii) constructing the sample cross-covariance matrix from the ensemble of model parameters and predicted responses. Once the cross-covariances between model parameters and predicted data are obtained, they are used to update model parameters from field data. The exact form of the update is obtained by linearly combining the initial (i.e., prior) ensemble of model parameters with the ensemble of (perturbed) observations, while imposing certain requirements on the solution (namely, unbiasedness and minimum estimation error variance). The success of EnKF in subsurface flow model calibration is largely due to the prevalent spatial continuity in rock property distribution and the strong correlations that exist between flow responses and formation hydraulic properties. In this section, we briefly review the EnKF formulation and highlight its main properties. A detailed treatment of the topic can be found in (Evensen et al., 2007), while a review of the application of the method to subsurface flow model calibration is given in (Aanonsen et al., 2009). The EnKF (Evensen, 1994) is a Monte-Carlo-based approximation of the Kalman filter (KF), which is a sequential Bayesian estimation method (Gelb, 1970). The KF is a linear state estimation method that provides optimal results (exact posterior distribution) for linear dynamical models in which the uncertain elements are described using multi-Gaussian statistics. The filter linearly combines prior estimates of the states with the likelihood function of the observed quantities to obtain the posterior state estimate. The weight of the linear combination, known as the Kalman gain, is the normalized (by predicted observation covariance) cross-covariance between the state variables and observations. The linearity assumption significantly simplifies error (covariance) propagation in the dynamical system while the Gaussianity assumption enables the use of mean and covariance to fully characterize the distribution. Developed to extend the use of KF to nonlinear problems, the EnKF provides a computationally efficient approximation for obtaining the mean and covariance needed in the regular linear KF update. It uses an ensemble representation of uncertain variables (by sampling from their PDF), which allows for nonlinear propagation of state variables (and possibly model parameters) forward in time. Using a similar update to KF, the ensemble version uses sample-approximated Kalman gain to estimate the states/parameters sequentially whenever measurements become available (Evensen, 1994; 2004). For completeness, we present an overview of the EnKF formulation next. We denote a general non-linear state-space model as:

(tk ) = fk ( (tk 1)) + wk d (tk ) = hk ( (tk )) + vk

T i ,

i

f

(tk ) = fk (

u i (tk 1))

(3)

where i f (tk ) denotes the ith predicted state realization at time tk and th u i (tk 1) refers to the i updated state realization at the previous time step. We note that model error (or process noise) is not included in the forecast step. One reason for excluding model errors in the formulation is that characterizing modeling errors is very complex. In addition, in subsurface flow model calibration literature, the uncertainty in describing model inputs is considered to be more dominant. Therefore, it is assumed that the errors related to accurate description of the physical processes are not as significant as those that result from uncertain input parameters. Excluding the model error results in compliance with the physical conservation principles that are used to derive the nonlinear state propagation model. However, it is important to note that several other inputs into the simulation model (e.g., relatively permeability functions, capillary pressure model, and other fluid properties) can be uncertain and involve potentially significant errors that may need to be accounted for. Since rock properties are assumed to be time-independent parameters, they remain the same during the forecast stage of the filter. When measurements become available, the EnKF update process is carried out to linearly combine forecast states with observations. The standard form of the EnKF update can be expressed as (Evensen et al., 2007):

x iu (tk ) = x i f (tk ) + Kk [dobs, i (tk )

Hx i f (tk )],

(4)

th

where dobs, i (tk ) represents the i realization of the perturbed measurements at time tk ; H is the measurement matrix with only zero and one elements as its entries (used to extract the measurements from the augmented state vector.), K is the Kalman gain matrix and is expressed as:

Kk = Ckf H T [HCkf H T + Rk ]

(5)

1

in which is the covariance matrix of the forecasting state computed as:

Ckf

Ckf =

Xkf (Xkf )T N 1

Xkf = Xkf

(1)

X¯ kf = Xkf

and is

(6)

Xkf 1N / N = Xkf (I

Xkf = {X1f (tk ), X2f (tk ), …, XNf (tk )}

where (tk ) is the system state vector at time tk , fk is a non-linear state propagation function with the additive model error , hk is a non-linear measurement operator with additive measurement error vk , and d (tk ) is the measurement vector at time tk . The EnKF filtering process is initialized with an ensemble of N initial (often augmented) state vectors = { 1, 2 , …, N } that are sampled from their respective probability density functions prior to observing any dynamic measurements. The ith augmented state vector, which is constructed by adding parameters and

f

n× N

1N / N )

(7) (8)

Here, n is the dimension of the augmented forecast state vector x i f , N ×N N × N is a matrix with all elements equal to 1, I 1N denotes the identity matrix, and Rk is the observation error covariance matrix that can be computed from the observation perturbation matrix Ek = { 1 (tk ), 2 (tk ), …, N (tk )} as follows

Rk =

201

Ek EkT N 1

(9)

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Fig. 1. Schematic representation of the EnKF forecast and update steps. The forecast step involves non-linear propagation of a finite number of state realizations to generate an ensemble of predicted observations; the analysis step updates each individual predicted state/parameter realization by linearly combing them with observed quantities using a Kalman gain matrix.

The perturbed measurements can be expressed as N dobs, i (tk ) = d (tk ) + i (tk ) , where d (tk ) is the measurement vector obtained at time tk . Various implementation of the EnKF update have been proposed in the literature. In this paper, we use the square root filter (SRF) that was proposed by (Evensen, 2004) to estimate the static parameters (i.e., porosity, horizontal and vertical permeability maps) of a CO2 storage aquifer model. The measurements used in our examples consist of the pressure and CO2 saturation at the existing wells, as well as the timelapse seismic data. The augmented state vector in our experiments can be written as:

{ , kh, k v , dBHP , d Sg , dseismic}

implementation reduces the dimension of the problem from the state size (typically in the order of ∼106–7) to the ensemble size (usually very small ∼102). As such, the updated parameters reside in a low-dimensional subspace defined by the prior parameter realizations. This property affects the EnKF robustness to bias and uncertainty in the prior model realizations (Jafarpour and McLaughlin, 2009; Jafarpour and Tarrahi, 2011). 3. Geologic model for CO2 storage We use a series of field-scale CO2 storage experiments to evaluate the EnKF performance and the relative significance of different data types. The structural model is borrowed from the PUNQ-S3 benchmark model and represents a folding geologic formation that forms an anticline (Floris et al., 2001). We have refined the grid system in the original PUNQ-S3 model to obtain a high-resolution 3D model to achieve a higher simulation accuracy. The heterogeneity in the properties of the storage formation is modelled using correlated Gaussian random fields (GRF). To generate correlated realizations of porosity, horizontal and vertical permeability fields, first a GRF conditioned on porosity data is generated based on a set of specified variogram model parameters, using the SGSIM algorithm (Deutsch and Journel, 1998). This model is then used as secondary data with a specified correlation coefficient (0.8 in this paper) between porosity and horizontal permeability to generate a horizontal permeability (conditioned on hard data) using sequential Gaussian co-simulation (Deutsch and Journel, 1998). The same procedure is repeated to obtain the corresponding conditional vertical permeability map (for which we also used correlation coefficient 0.8). The reference model for porosity and permeability are also generated using the sequential Gaussian simulation (SGSIM) algorithm (Deutsch and Journel, 1998). The variogram model parameters and the statistics used in generating the reference model are listed in Tables 1a and 1b. In the reference model, layers 1–3, 7–9 and 13–15 have similar characteristics. However, layers 4–6 have much lower hydraulic properties, while, layers 10–12 have intermediate properties. Fig. 2(a) shows the injection and monitoring well locations and the map of gas saturation after 4 years. Fig. 2(b) displays, from left to right, the horizontal and vertical permeability, and the porosity for the reference model. We note that in this study we treat variogram model parameters (anisotropy direction and ratio) as random variables to reflect the uncertainty in the initial ensemble of formation properties. The geologic CO2 storage simulations presented in this paper are carried out using the CO2STORE option of the Eclipse ® subsurface simulation package (E300 module) (ECLIPSE, 2012). The model assumes three phases (CO2-rich gas phase, water-rich liquid phase, and solid phase representing salt precipitations) with water, gas, and salt (NACL,

(10)

where is porosity, kh and k v stand for the horizontal and vertical permeability, respectively, dBHP and d Sg refer to the predicted pressure and CO2 saturation measurements at the monitoring wells, and dseismic refers to the predicted 4D seismic data. With this definition of the state vector, it is possible to update the uncertain static parameters (kh and k v ) using pressure, CO2 saturation and time-lapse seismic data. The EnKF update with high-dimensional observations can suffer from ensemble collapse and filter divergence issues if the update does not account for the true rank of the variables involved (Skjervheim et al., 2007). To integrate high-dimensional and correlated time-lapse saturation data, we follow the approach proposed in (Skjervheim et al., 2007; Tarrahi et al., 2015), where low-rank representations of the data is used for assimilation. The subspace used for this representation is defined by the left singular vectors of the predicted observations matrix (for details see (Skjervheim et al., 2007; Tarrahi et al., 2015)). Denoting the perturbed predicted observation matrix as and collecting its leading left singular vectors into U = [u1 u2 … uN ], the projected perturbed observation ensemble Dp is calculated as Dp = U TD . In our experiments, to simulate the errors involved in inverting the 4D seismic data to timelapse saturation maps, we have added uncorrelated (white) noise to the projection coefficients of the low-rank representation of the simulated time-lapse saturation observations. This approach is equivalent to adding correlated (red) noise to the spatial representation of the timelapse saturation observations. In our EnKF implementation, the updated states of the system are obtained by first estimating aquifer properties, and then using the updated properties in reservoir simulation to derive the updated states (from the very beginning of the simulation), a process known as the restart option (Wen and Chen, 2006; Hendricks Franssen and Kinzelbach, 2008; González-Nicolás et al., 2015). Evensen (2004) shows that the EnKF update equation for each state/parameter realization can be reformulated and expressed as an expansion over the predicted ensemble members. This implies that the filter 202

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

from the distribution specified in Table 3. The resulting realization is used in the SGSIM algorithm to generate initial realizations of porosity and permeabilities. This process is repeated 100 times, each time with difference variogram model parameters, to generate 100 initial realizations for the EnKF implementation. Fig. 3 shows two sample realizations, the mean and variance of the initial ensemble of kh , k v and . The measurements from the wells were corrupted by adding an uncorrelated Gaussian noise with a standard deviation of 20% of each measurement. For integration of 4D seismic, we added correlated noise to vertically-averaged saturation difference maps. In general, the relationship between saturation and time-lapse seismic (amplitude differences) is complex and the quality of the inverted saturation results is field-specific. The seismic response to an initial increase in CO2 saturation (0–10% saturation) is known to be very strong. However, the sensitivity diminishes rapidly as CO2 saturation increases beyond this range (Roach et al., 2017). Furthermore, the vertical resolution of seismic data depends on the thickness and depth of the formation (Roach et al., 2017), with a better chance of detection with vertical resolution for shallower and thicker reservoirs. Another important factor is the contrast between CO2 and brine properties, which results in a strong seismic response as CO2 substitutes brine in the storage formation. To account for low-resolution of seismic data and to mimic realistic settings, instead of using saturation differences for each layer, we represent time-lapse seismic data with vertically-averaged CO2 saturation differences and add correlated noise to account for the uncertainty in processing and interpretation of time-lapse seismic data. This averaging approach tends to be a conservative assumption as the loss of vertical resolution in time-lapse seismic data reduces its information content and assimilation impact. We note that in some cases, it may be more appropriate to assume that the time-lapse seismic only indicates the evolution of the CO2 plume within the top layer of the storage formation, i.e., beneath the cap rock (Cavanagh, 2013; Chadwick and Noy, 2010; Zhu et al., 2015). However, in relatively thick reservoirs, it is possible to detect CO2 saturation below the top layers, which is why we have used vertically averaged saturation differences in this work. The summary of data assimilation information is provided in Table 4.

Table 1a Information about the porosity and permeability of the Reference model. Statistical Information Layers Porosity 1–3 4–6 7–9 10–12 13–15 log10 kh 1–3 4–6 7–9 10–12 13–15 log10 kv 1–3 4–6 7–9 10–12 13–15

Avg.

Std.

Maxup.

Maxlow.

Limup.

Limlow.

0.14 0.08 0.14 0.1 0.14

0.11 0.04 0.11 0.06 0.11

0.3 0.17 0.3 0.22 0.3

0.01 0.01 0.01 0.01 0.01

0.25 0.12 0.25 0.17 0.25

0.05 0.03 0.05 0.04 0.05

2 1.4 2 1.6 2

1 0.6 1 0.8 1

3 2.3 3 2.7 3

−1 −1 −1 −1 −1

2.7 2 2.7 2.5 2.7

0 0 0 0 0

1.6 0.8 1.6 1.1 1.6

1 0.6 1 0.8 1

2.7 1.7 2.7 2 2.7

−1 −1 −1 −1 −1

2.3 1.4 2.3 1.7 2.3

0 0 0 0 0

Table 1b Variogram used to generate the Reference model. Variogram Model Parameters Layers

Correlength Length

Angle

Anisotropy

1–3 4–6 7–9 10–12 13–15

60 30 60 30 60

−15 0 −15 0 −15

3 1 3 1 3

CACL2, and CACO3) components. In the model, the partitioning of the gas and liquid phases follows the procedure proposed by Spycher and Pruess (2005). The detailed description of the simulation parameters in our experiments is provided in Tables 2a and 2b. Two injection wells and six monitoring wells are placed in the formation (Fig. 2(a)). The injectors are perforated only in the bottom layer (15th layer). Approximately 1.2 pore volume (in standard condition) of CO2 is injected into the formation during 10 years of injection. The total simulation time is 100 years, with 90 years of post-injection period to investigate how CO2 plume migrates in the formation. We consider two years of data assimilation and 98 years of forecast for performance evaluation. For data assimilation, we use four sets of well pressure and Sg measurements (collected at six months intervals) and a single time-lapse seismic data that represents the difference in CO2 saturation at the initial condition (base survey) and after two years (monitor survey). The update with time-lapse data is performed only once after two years when the monitor survey is obtained. In practice, the frequency of timelapse seismic data is limited by the level of detectable changes in the formation and the high cost of data acquisition.

4.1. Sensitivity analysis Prior to data assimilation, we examine the sensitivity of model responses (saturation and pressure) to static model parameters ( , kh and k v ) that will be estimated. An exhaustive sensitivity analysis that accounts for the coupling among parameters is computationally demanding to perform. We consider a few relevant cases in our sensitivity analysis to gain insight about the information in different parameter types. First, we consider a homogeneous base case, where , kh and k v are constant in all grid cells. These constant values are obtained by taking the field average of the respective parameters in the reference model. For each parameter, flow simulations are run for two cases, each generated by adding ± 5% perturbation to the base case to approximate the sensitivity of the model responses at different time steps. Similar experiments are performed for kh and k v , and the results are listed in the first row of Table 5, which show that both saturation and pressure data exhibit slightly more sensitivity to the porosity and horizontal permeability. We also considered another set of experiments to investigate the sensitivity information for a layered model, where the base case has constant model parameters in each layer. Following a similar procedure, the sensitivity analysis results (shown in the second row of Table 5) suggest the model responses, including saturation and pressure, are more sensitive to porosity and horizontal permeability. Several authors have studied the importance of aquifer heterogeneity on CO2 injectivity, displacement, and storage efficiency of the injected CO2 (e.g., (Kumar et al., 2005; Jahangiri and Zhang, 2011; Han et al., 2014; Tian et al., 2014, 2016)). We also performed a similar sensitivity

4. Experiments and results We consider the estimation of log10 (kh) , log10 (k v ) and as uncertain hydraulic properties of the aquifer using the square-root filter implementation in (Evensen, 2004). Notice that in real applications, the true variogram model and it parameters are usually unknown. In this study, we assume the parameters of the variogram models for generating the initial ensemble of porosity and permeabilities are uncertain and use random variables with the distributions listed in Table 3 to describe them. For each parameter combination (correlation length, azimuth angle, and anisotropy ratio), a random sample is generated 203

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Fig. 2. Reference model and well configuration: (a) well locations and CO2 saturation distribution after four years of injection; (b) from left to right, horizontal and vertical permeability (kh, k v ) and porosity ( ) of the Reference model. Table 2a Summary of simulation model and fluid properties.

Table 2b CO2 injection rates.

Parameter

Value

Time

INJ-1 (105m3/ day )

INJ-2 (105m3/ day )

Model Size (Number of Grids) Grid Dimenssion (X, Y and Z Direction) Phases Reservoir Temperature Datum Pressure Datum Depth Rock Compressibility

57 × 84 × 15 60 m × 60 m × 2 m Water/Gas 34 234 Bar 2355 m

0–6 months 7–12 months 13–18 months 19–24 months

3.2002 2.1817 2.7657 2.7176

2.3998 3.4183 2.8343 2.8824

Reference Density of Solid Phase Components Number of Injectors Total Injection Time Total Injection Volume Injection Well Constraints Number of Monitoring Wells Measurements Total Simulation Time

5 × 10 5 PV/PV/psi−1 2650 kg/m3 H2 O / CO2 / NaCl 2 10 years 1.2 Pore Volume Injection Rates 6 BHP / Sg / 4D Seismic 100 years

Table 3 Uncertain variogram model parameters (Gaussian distribution).

analysis with respect to the heterogeneity in the reference model. In this case, we uniformly increased and decreased the property values by applying a 5% multiplier to all grids and used the results to compute the sensitivities. The results are listed in the third row of Table 5 and show similar outcomes. Similar experiments with heterogeneity borrowed from other realizations provided additional support for these

Layers

Correlength Length

Angle

Anisotropy

Avg. Std. Min. Max.

45 15 5 90

−15 15 −60 30

2 0.3 1 3

observations (Table 5). The perturbation-based sensitivity analysis in Table 5 represents local sensitivity information and the relative importance of each data in estimating different parameters. We also computed the sensitivity information based on the correlation between model responses and parameters in the heterogeneous case (results not shown). 204

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Fig. 3. Summary representation of the initial ensemble of uncertain hydraulic aquifer properties, kh (left), k v (middle), and (right). Two initial samples, as well as the mean and variance maps are shown in Rows 1 through 4, respectively. The realizations are conditioned on static well data.

In Eqs. (4) and (5), model parameters are updated by using their cross-covariance with the model response. In our study, the generated initial ensemble of permeability and porosity maps have a correlation coefficient of 0.8. This implies that if model responses are correlated with one parameter (e.g., porosity), they are also correlated with the other two parameters. Consequently, the updates applied to model parameters tend to be similar. The model responses, including the saturation and pressure, to the initial ensemble of parameters are shown in Fig. 4. The total variance map of the predicted Sg in each layer, over 100 years of simulation, is calculated by taking the sum of variance values in all grid cells, layer by layer. The red solid line represents the time when CO2 injection stops

Table 4 Data assimilation information. Updated Parameters

log10 (kh) , log10 (kv ) and

Number of Updates

4 Sg /BHP and 1 Seismic

Update Frequency

Subspace Representation of Seismic Data EnKF Update Method Number of Measurements (BHP, Sg ) at Each Update

Every 180 days for BHP and Sg Every 720 days for Seismic Dim = 80 Square Root Filter (SRF) 90 (6 wells in 15 layers)

205

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Table 5 Sensitivity analysis results. Pressure

Sg

Homogeneous Constant Layers Reference Model

log10 (kh)

log10 (kv )

1.6282 1.8104 0.2845

0.2793 0.8707 0.3175

2.0088 3.225 2.1185

log10 (kh)

log10 (kv )

123.5406 51.0599 88.1593

46.171 18.494 7.3937

174.4888 243.4025 155.233

4.2. Data assimilation with EnKF We present three experiments (Table 5) to investigate the role of different measurement types in both parameter estimation and future performance prediction. In all cases, log10 (kh) , log10 (k v ) and are updated as uncertain parameters, represented with an ensemble size of 100. In Experiment 1, only well pressure data are integrated. In Experiment 2, well pressure and CO2 saturation measurements are assimilated, while in Experiment 3 well pressure, CO2 saturation and time-lapse seismic data are used to update model parameters. To illustrate the parameter distributions before and after updates, we show areal and cross-sectional views of the property maps. Fig. 5 shows the two cross-sections that are selected for showing the property distributions (each passing through one injection well). Since the estimation results for log10 (k v ) are similar to those of log10 (kh) , we only show the final (updated) ensemble of log10 (kh) and (Figs. 6–9, respectively). In each figure, the first row shows the reference map for the parameter of interest, while the second to fifth rows depict two sample realizations, the mean and variance maps corresponding to the initial ensemble (first column) and the final ensembles for Experiments 1 to 3 (second to fourth columns, respectively). In Experiments 1 and 2, the updated mean and variance maps show minimal changes when compared with their initial ensemble. In Experiment 3, however, the updated mean and variance maps show more significant changes. Overall, the updated mean from Experiment 3 shows significant improvement in permeability distribution. The results in Figs. 6–9 show only very small uncertainty reduction in the final ensemble of porosity and permeability in Experiment 1 and 2. The small updates (and variance reduction) in Experiments 1 and 2 indicate limited constraining power for well pressure and saturation data. We also note that the CO2 saturation does not reach some of the monitoring wells during the first two years (data assimilation period). This observation points to two important considerations regarding well data for monitoring: first, it is very important where the wells are placed to maximize their information content; second, data from wells at fixed locations are not very effective for monitoring the CO2 plume evolution as they cannot track the evolution of the CO2 front. The performance of data assimilation for each experiment is quantified and listed in Table 6, where the reductions (relative to the initial ensemble) in total variance of model parameters after the last update in each experiment are reported. The total variance is obtained by summing the variance in all grid cells. The results from Experiment 3, in which 4D seismic data is integrated, show larger variance reduction. Table 7 lists the relative reduction in total RMSE of model parameters after data assimilation in Experiments 1–3. The total RMSE reduction in Experiment 3 is also significantly larger than the RMSE reductions in Experiments 1 and 2, implying more significant error reduction for updated parameters in Experiment 3. These observations confirm that, compared with the pressure and saturation data in scattered wells, time-lapse seismic data can play a more significant role in constraining model parameters. The relative reduction in the total RMSE (Fig. 10, top) and total variance (Fig. 10, bottom) of the predicted Sg for the final ensembles in Experiments 1–3 show the superior results obtained in Experiment 3. In

Fig. 4. Total variance of CO2 saturation (Sg ) predictions in each layer at different time steps. The red solid line refers to Year 10, when CO2 injection stops. Note: after Year 10 the time-axis is not linear. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Fig. 5. Top:two cross-sectional slices (along x-axis) in the model that pass through the injection wells; bottom: areal view of the two cross-sections.

(10 years). Note that the x-axis for Years 10–100 in Fig. 4 is not plotted in linear scale. It can be observed that during the injection of CO2, there is little difference in the total variance of the layers at any given time, especially in the first 5 years. After the injection period, however, there is a clear separation in the total variance. The response of the upper layers in the formation is very different from (i.e., higher than) the response in lower layers. This can be explained by the change in the driving forces in migration of the CO2 plume. During injection, pressure gradient and buoyancy forces both play important roles in displacing the CO2 plume. Once the injection stops, however, the migration of the CO2 plume is dominated by the gravity force. Consequently, the CO2 plume movement is primarily in upward direction. Hence, the major changes in the aquifer occur in the upper layers of the formation with the lower parts containing small (residual) amounts of CO2. The change in driving forces after injection also implies that the prediction of Sg in the lower part of the formation is less sensitive to rock properties in those layers. Hence, updates applied to the properties in the lower regions have smaller impacts on the predicted Sg .

206

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Fig. 6. The initial and final (updated) ensemble of log10 kh in Experiments 1–3.

Fig. 10, after Year 10 (when CO2 injection stops) a sudden change is observed in the behavior of these plots. To explain this behavior, Fig. 11 shows the total RMSE and total variance of Sg predictions in each layer and at all time steps (for 100 years; the red line indicates Year 10 when CO2 injection stops). As can be seen, due to predominantly buoyancy forces after Year 10, higher values of RMSE and variance are concentrated in upper layers. Referring to Fig. 7 (and Fig. 9), the crosssectional view of the updated ensemble of model parameters indicates significant variance reduction (and improved parameter estimates) in Experiment 3. The improvements in the estimated properties for top layers exhibit themselves as reduced RMSE and variance in long-term saturation predictions. The predicted CO2 storage (in residual and dissolved forms) with the initial and update model parameters from different experiments are summarized in Figs. 12 and 13. Each boxplot in Fig. 12 shows, the interquartile ranges (q3 q1), where q3 and q1 are the 75th and 25th percentiles, respectively, for the amount of CO2 storage after years 20 to 100. The dash inside each rectangle shows the median and the whiskers

above and below the box show the location of q1 1.5(q3 q1) and q3 + 1.5(q3 q1) , respectively. For normally distributed variables, the interval between the two whiskers covers approximately 99.3% of the values. The red stars refer to the values from the reference model while the small black stars show the mean values of each ensemble predictions. The predicted storage from Experiment 3 show pronounced improvements, both in terms of accuracy and confidence level. Fig. 13 shows the reduction (relative to the initial ensemble) in RMSE for the predicted CO2 storage in Experiments 1–3. The final ensemble from Experiment 3 also has the largest relative RMSE reduction (more than 40%) in predicting the residual and dissolved CO2, confirming the observations made about the improvement in estimated permeability and porosity maps. 4.3. Formation with natural fractures A major risk in GCS projects is the CO2 leakage from the formation through existing or created pathways (wells, faults, fractures and 207

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Fig. 7. Two cross-sectional views (see Fig. 5) of the initial and final ensemble of log10 kh in Experiments 1–3.

compromised caprocks). Formations with natural fractures present an increased risk of leakage, especially when such fractures exist in the caprock or can propagate toward it. While faults and large fractures may be detected in 3D seismic images, the presence of natural fractures can only be estimated from well logs, outcrop, or core samples. Another possibility is creation of induced fractures in response to injection-related stress variations. An important question is whether assimilation of monitoring data can provide insight about the presence of undesirable flow features or pathways (such as fractures) in the formation. In this section, we investigate the possibility of detecting the presence of unanticipated flow features, such as natural fractures, from assimilating CO2 storage monitoring data. We present two experiments: Experiment A and B. In Experiment A the observed data is from a subsurface model with no fractures while in Experiment B, the data is from a reservoir with natural fractures. In both experiments, the initial ensemble of model parameters is the same and does not include any fractures. All parameters used in the EnKF implementation remain identical. Fig. 14 shows the permeability and porosity of the reference model with fracture networks (dark colors) in Experiment B. The reference model in Experiment A has the same property distributions as in Experiment B, but without the natural fractures. The use of EnKF in these experiments is not intended to recover the fracture distributions, but mainly to investigate how the filter response changes due to presence of fracturs. The results of EnKF for these experiments are shown in Fig. 15. In this figure, the top row shows the reference log10 (kh) model and the

following rows show two samples, the mean and variance maps of the initial ensemble (left) and final ensembles in Experiment A (middle) and Experiment B (right). While in Experiment A there are no indications of the existence of extreme flow features in the permeability maps, the results from the final ensemble in Experiment B show areas with extremely high and out-of-range permeability values, pointing to an unusual filter behavior. While in our experiment this behavior is known to be caused by the presence of fractures, in practical applications such a response from a filter typically prompts modelers about a systematic issue in the design and implementation of the filter. To identify whether this behavior can be attributed to possible inadequacy of the initial ensemble, one can run a synthetic test case where the true model is taken to be one of the initial realizations. If application of the filter to this synthetic test case does not yield similar unexpected outcomes, the issue can be linked to inconsistency between prior ensemble and the reference model. We note that direct application of EnKF for estimation of discrete fractures is not advisable even when the initial ensemble contains such features. The reason for this is that EnKF is a continuous estimation method that has optimal performance in problems involving multiGaussian statistics, whereas fractures represent discrete and highly nonGaussian features. Compared with the reference model, the high permeability zones that are identified in Experiment B are mainly in areas where natural fracture networks are present in the reference model. The results from the two experiments suggest that EnKF may hint to the absence of important (extreme) features in the initial ensemble of 208

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Fig. 8. The initial and final (updated) ensemble of

model parameters if the observed data contains signatures of such anomalies. This provides an important feedback from data integration to correct the initial ensemble and possibly revise the inversion method if appropriate. While several approaches have been proposed in the literature to improve the EnKF performance for discrete non-Gaussian problems, treatment of discrete fractures is especially challenging due to the their complex geometric and flow attributes. In addition to data assimilation considerations, proper treatment of fractured models may require introducing additional complexities to the forward flow simulation, including non-Darcy flow effects, geomechanical considerations, and stress-dependent changes in fracture hydraulic properties.

in Experiments 1–3.

results suggest that without including time-lapse seismic data, it is less likely to obtain additional insight about the evolution of the CO2 plume in the aquifer as model calibration against limited well data may result in minor changes/improvements. In our experiments, data from five monitoring and two injection wells were used. In practice, it is highly unlikely to have many monitoring wells. Furthermore, monitoring the evolution of the CO2 plume with well locations that are fixed in space is not effective even when several wells are available. Since limited number of wells are expected in geologic CO2 storage, it is important to maximize the utility of those wells by collecting more informative monitoring information, e.g., through cross-well imaging and monitoring techniques. Whenever possible, time-lapse surface seismic data should be combined with cross-well seismic (Ajo-Franklin et al., 2013; Butsch et al., 2013) and other monitoring information to enhance the constraining power of the resulting data. Other types of monitoring data, such as passive microseismic monitoring (Bauer et al., 2014; Finley, 2014), gravimetry measurement (Alnes et al., 2011), deformation and ground-level elevation changes using Interferometric synthetic

5. Discussion While the presented formulation and results in this paper are based on simplifying assumptions, they offer a starting point to consider EnKF for monitoring data assimilation in geologic CO2 application. The 209

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Fig. 9. Two cross-sectional views (see Fig. 5) of the initial and final (updated) ensemble of

effects. Coupled modeling of thermal-hydrological-mechanical-chemical (THMC) processes and dynamic changes in rock/fracture properties during and after injection are necessary for building the correlations among multi-modal data and physical rock properties. While preliminary results from applying the EnKF with coupled flow and geomechanics models to integrate microseismic monitoring data in enhanced geothermal systems (EGS) have been promising (Tarrahi et al., 2015), it is important to note that multi-physics models introduce additional complexity and nonlinearity that may have an impact on the performance of EnKF. In this study, to represent time-lapse seismic instead of using seismic amplitude changes we used noise-corrupted vertically-averaged saturation differences in time. This was done to avoid seismic modeling and focus primarily on the data assimilation with EnKF. Addition of noise was used to account for the errors that may be introduced during seismic processing. While it is common in reservoir engineering to use saturation differences to represent time-lapse seismic information, a more comprehensive approach is to combine flow simulation with seismic wave propagation and rock physics models to assimilate timelapse seismic amplitude differences. Another important assumption that we used in representing seismic data was vertical averaging for the saturation differences. This choice was motivated by the limited vertical resolution of the time-lapse seismic data. It is important to emphasize that detection and representation of CO2 saturation with time-lapse seismic is complex and involves considerable uncertainty. In general, the accuracy and repeatability of data, the geometry of the CO2

Table 6 Reduction in total variance after data assimilation in Experiments 1–3.

Final 1 (Update 4) Final 2 (Update 4) Final 3 (Update 4)

kh

kv

3.95% 4.95% 39.30%

8.40% 9.73% 41.64%

2.39% 3.50% 39.02%

Table 7 Reduction in total RMSE after data assimilation in Experiments 1–3.

Final 1 (Update 4) Final 2 (Update 4) Final 3 (Update 4)

kh

kv

0.55% 1.12% 9.45%

2.00% 2.95% 10.64%

in Experiments 1–3.

0.32% 0.55% 8.91%

aperture radar (InSAR) data (Ringrose et al., 2013; White et al., 2014), as well as near-well geophysical monitoring (Sakurai et al., 2005; Hovorka et al., 2006; Sato et al., 2011) can also provide complementary information that can be integrated into reservoir models. While EnKF has shown great promise in handling multiple data types and simultaneously estimating several types of parameters, a main difficulty in integration of multi-modal data is the proper representation and simulation of the related multi-physics processes, and their coupling 210

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

reported. A major difficulty is related to the degradation of imaging quality in deeper layers, which limits the constraining power of timelapse seismic. This issue has also prompted recent studies to focus on imaging and characterizing the topmost layer of CO2 beneath the formation top seal (Cavanagh, 2013; Chadwick and Noy, 2010; Zhu et al., 2015). This approach is motivated by the observation that the longterm evolution of the CO2 plume is controlled by the top layer of the formation as a large portion of the injected CO2 will rise to top layer and spread laterally. However, the buoyancy behavior of the plume may also be influenced by reservoir properties and field-specific heterogeneity, which can introduce significant uncertainty into the interpretation process. Another important consideration in representing time-lapse seismic data is data dimensionality. In general, the most informative part of the seismic data is related to CO2 front-tracking. This implies that timedifferential seismic images are likely to have overwhelmingly zero/insignificant components (areas with no or small changes). Therefore, the resulting time-lapse seismic images are highly redundant and can be represented in low-dimensional subspaces. This is an important issue in implementing the EnKF as without a low-rank representation of timelapse seismic the filter performance can be significantly degraded. However, the choice of a low-dimensional space to effectively represent time-lapse seismic data is an open area for investigation. While we used the leading left-singular vectors of the perturbed predicted (verticallyaveraged) saturation differences to construct a low-rank data representation prior to assimilation, other descriptions that focus on CO2 front-tracking may also be used (Leeuwenburgh and Arts, 2014; Trani et al., 2017). Lastly, the use of uncertain variogram model parameters in our study was motivated by the limited state of knowledge about these parameters in a saline aquifer. During model calibration, we only estimated the spatial distribution of rock properties in the reservoir, without updating geostatistical (e.g., variogram) model parameters. It may be possible, however, to gain insight into the geostatistical

Fig. 10. Relative reduction (compared with the initial ensemble) in total RMSE (top) and variance (bottom) of the predicted Sg in Experiments 1–3.

accumulation, reflectivity and properties of the CO2 itself, as well as reservoir thickness and depth play important roles in detecting changes in CO2 saturation. While efforts to combine flow simulation with timelapse seismic monitoring data have resulted in important insight about the evolution of the CO2 plume (Lindeberg et al., 2001; Van der Meer et al., 2001), significant levels of uncertainty and non-uniqueness are

Fig. 11. Total variance (left column) and RMSE (right column) of CO2 saturation (Sg ) predictions in each layer at different times. The red line shows Year 10, when CO2 injection stops. Note: after Year 10 the time-axis is not linear. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). 211

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Fig. 12. Prediction of stored CO2 (through residual and dissolution trapping) using the initial ensemble and the three final (updated) ensembles from Experiment 1–3.

parameters of the reservoir properties by jointly inverting for those parameters as well as the spatial distribution of the formation hydraulic properties. Methods for joint estimation of rock properties and the statistical parameters of their distribution have been reported in the literature (Rubin et al., 2010; Jafarpour and Tarrahi, 2011; Jardani et al., 2012; Hansen et al., 2012; Laloy et al., 2015). While some studies have found that a combination of data scarcity and non-unique and complex relation between flow data and geostatistical parameters (e.g. variogram sill and range) limit the ability to constrain these parameters during flow model calibration (Jafarpour and Tarrahi, 2011), others have reported successful inference of geostatistical parameters from flow data for simple 2D examples (Rubin et al., 2010; Jardani et al., 2012; Hansen et al., 2012; Laloy et al., 2015). An alternative method for model calibration under uncertainty in prior geostatistical scenarios is to initially include multiple plausible geologic scenario and use dynamic data to first identify consistent prior geologic scenario(s) (Khaninezhad and Jafarpour, 2014; Golmohammadi and Jafarpour, 2016).

knowledge about the properties of saline aquifers in realistic situations, the initial ensemble of formation hydraulic properties was obtained by assuming significant uncertainty about the correlation structure and anisotropy direction and length, (i.e. parameters of the variogram model). Despite the low-resolution and information content of timelapse saturations (due to vertical averaging), and highly uncertain nature of the initial ensemble used, the filter was able to estimate the general trend in the permeability and porosity distribution, especially in areas traversed by the CO2 plume. Several interesting observations were made in our experiments. First, the spatial and temporal distribution of the sensitivity information exhibited different behavior during versus after CO2 injection. During injection, where both pressure gradient (together with permeability heterogeneity) and buoyancy forces are important, the CO2 plume migrates both laterally and vertically in the aquifer. As a result, the plume migration (and CO2 saturation) shows sensitivity to rock formation properties in all layers. However, after injection stops, the buoyancy becomes the dominant driving force in displacing the CO2 plume, with heterogeneity and flow barriers being the main mechanisms for lateral movement of the plume. In this case, the plume migration becomes more sensitive to formation properties mainly in the upper layers. This difference in sensitivity information also affects the EnKF update and the resulting long-term prediction of the CO2 plume. While poor estimates of rock properties in top layers may not be detrimental in initial stages, they become important at later times. Another observation is related to well data conditioning and its impact on initial displacement of the CO2 plume. Our results show that conditioning the initial ensemble of rock permeability and porosity models on hard data from injection wells leads to similar CO2 saturation distributions around these wells, especially at earlier stages. As CO2 moves laterally in the formation, the ensemble spread for saturation distribution tends to increase with growing variation in rock properties. We also studied the effect of inconsistency in the prior model. Specifically, when the initial ensemble of aquifer properties does not include important flow features, such as natural fractures. The results from our experiments suggest that when fracture networks exist but are

6. Conclusion In this paper, we investigated the feasibility of applying the EnKF as a data assimilation approach to characterize geologic CO2 storage aquifers from monitoring measurements, including scattered monitoring well pressure and saturation data, and time-lapse seismic data. Time-lapse saturation maps over a two-year interval were vertically averaged across all layers of the simulation model, and corrupted with spatially correlated noise, to represent 4D seismic data. A perturbationbased sensitivity analysis was used to study the effect of the aquifer porosity, horizontal and vertical permeability distributions on model responses (including pressure and CO2 saturation). The results in our experiments show that CO2 saturation and pressure distributions are sensitive to both porosity and permeability distributions. The EnKF data assimilation method was applied to a refined version of the PUNQ-S3 model, which represents a synthetic 3D case study based on an anticline geologic structure. To represent data scarcity and lack of reliable 212

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

Fig. 13. Plot of relative reduction (compared with the initial ensemble) in the RMSE of the predicted stored CO2 from the three final ensembles in Experiment 1–3. From top to bottom, residual, dissolved and total (residual + dissolved) storage are shown.

Fig. 15. The log10 (kh ) distribution for the initial ensemble (left), final ensembles in Exp. A (middle) and Exp. B. (right).

not included in the initial ensemble, the filter provides updates that are not consistent with the initial representation of rock properties, signaling a modeling issue that needs to be addressed. The use of fractures in these examples was motivated by the possibility of having natural fractures or injection-induced fractures in the storage formation. However, the presence of fractures in the aquifer may trigger complex issues to be considered in both forward modeling (e.g., flow and geomechanical effects) and inverse modeling (e.g., strong nonlinearity and

non-Gaussianity effects). Finally, in this paper we included time-lapse seismic data as time-differential saturation maps, without modeling/ processing the raw seismic data. To represent the low resolution of seismic data and the errors that are introduced in seismic data processing, we used vertically-averaged saturation maps and added spatially correlated noise to the saturation differences, respectively. These

Fig. 14. Reference model of horizontal and vertical permeability (log10 (kh ) , log10 (k v ) ), and porosity ( )in Exp. B with fracture networks (shown in dark grey). 213

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al.

assumptions were used to closely represent the information content of time-lapse seismic data without the complexities associated with rock physics and seismic models, which enabled us to focus primarily on the feasibility of using EnKF for data assimilation. Field application of geologic CO2 storage is expected to involve several complexities and challenges related to modeling, monitoring, and characterization of these systems, including those discussed above, that must be investigated and addressed to ensure successful implementation.

Finley, R.J., 2014. An overview of the Illinois Basin? Decatur project. Greenhouse GasesSci. Technol. 4, 571–579. Floris, F.J.T., Bush, M.D., Cuypers, M., Roggero, F., Syversveen, A.-R., 2001. Methods for quantifying the uncertainty of production forecasts. Pet. Geosci. 7, 87–96. Golmohammadi, A., Jafarpour, B., 2016. Simultaneous geologic scenario identification and flow model calibration with group-sparsity formulations. Adv. Water Resour. 92, 208–227. https://doi.org/10.1016/j.advwatres.2016.04.007. ISSN 0309-1708. González-Nicolás, A., Baù, D., Alzraiee, A., 2015. Detection of potential leakage pathways from geological carbon storage by fluid pressure data assimilation. Adv. Water Resour. 86, 366–384. Han, W.S., Kim, K.Y., Choung, S., Jeong, J., Jung, N.H., Park, E., 2014. Non-parametric simulations-based conditional stochastic predictions of geologic heterogeneities and leakage potentials for hypothetical CO2 sequestration sites. Environ. Earth Sci. 71 (6), 2739–2752. Hansen, T.M., Cordua, K.C., Mosegaard, K., 2012. Inverse problems with non-trivial priors: efficient solution through sequential Gibbs sampling. Comput. Geosci. 16, 593–611. https://doi.org/10.1007/s10596-011-9271-1. Hassanzadeh, H., Pooladi-Darvish, M., Keith, D.W., 2009. Accelerating CO2 dissolution in saline aquifers for geological storage-mechanistic and sensitivity studies. Energy Fuels 23, 3328–3336. Hattenbach, R., Wilson, M., Brown, K., 1998. Capture of carbon dioxide from coal combustion and its utilization for enhanced oil recovery. In: Eliasson, B., Riemer, P., Wokaun, A. (Eds.), Greenhouse Gas Control Technologies. August 30-October 2, 1998. 4th International Conference on Greenhouse Gas Control Technologies (GHGT4), Interlaken, Switzerland. Hendricks Franssen, H.J., Kinzelbach, W., 2008. Real-time groundwater flow modeling with the ensemble Kalman filter: joint estimation of states and parameters and the filter inbreeding problem. Water Resour. Res. 44, W09408. https://doi.org/10.1029/ 2007WR006505. Hovorka, S.D., Benson, S.M., Doughty, C., Freifeld, B.M., Sakurai, S., Daley, T.M., Kharaka, Y.K., Holtz, M.H., Trautz, R.C., Nance, H.S., Myer, L.R., Knauss, K.G., 2006. Measuring permanence of CO2 storage in saline formations: the Frio experiment. Environ. Geosci. 13, 105–121. Jacquard, P., 1965. Permeability distribution from field pressure data. Old Spe J. 5 (04), 281–294. Jafarpour, B., McLaughlin, D.B., 2009. Reservoir characterization with the discrete cosine transform. Spe J. 14 (01), 182–201. Jafarpour, B., Tarrahi, M., 2011. Assessing the performance of the ensemble Kalman filter for subsurface flow data integration under variogram uncertainty. Water Resour. Res. 47https://doi.org/10.1029/2010WR009090. W05537. Jahangiri, H.R., Zhang, D., 2011. Effect of spatial heterogeneity on plume distribution and dilution during CO 2 sequestration. Int. J. Greenhouse Gas Control. 5 (2), 281–293. Jahns, H.O., 1966. A rapid method for obtaining a two-dimensional reservoir description from well pressure response data. Old Spe J. 6 (04), 315–327. Jardani, A., Dupont, J.P., Revil, A., Massei, M., Fournier, M., Laignel, B., 2012. Geostatistical inverse modeling of the transmissivity field of a heterogeneous alluvial aquifer under tidal influence. J. Hydrol. 472, 287–300. Jenkins, C., Chadwick, A., Hovorka, S.D., 2015. The state of the art in monitoring and verification–Ten years on. Int. J. Greenhouse Gas Control 40, 312–349. https://doi. org/10.1016/j.ijggc.2015.05.009. Jessen, K., Sam-Olibale, L., Kovscek, A., Orr, F.M., 2001. Increasing CO2 Storage in Oil Recovery. First National Conference on Carbon Sequestration May 14-17, 2001. Juanes, R., Spiteri, E.J., Orr, F.M., Blunt, M.J., 2006. Impact of relative permeability hysteresis on geological CO2 storage. Water Resour. Res. (12), 42. Khaninezhad, M.R., Jafarpour, B., 2014. Prior model identification during subsurface flow data integration with adaptive sparse representation techniques. Comput. Geosci. 18 (1), 3–16. Kumar, A., Noh, M.H., Ozah, R.C., Pope, G.A., Bryant, S.L., Sepehrnoori, K., Lake, L.W., 2005. Reservoir Simulation of CO2 Storage in Aquifers. Society of Petroleum Engineershttps://doi.org/10.2118/89343-PA. Laloy, E., Linde, N., Jacques, D., Vrugt, J.A., 2015. Probabilistic inference of multiGaussian fields from indirect hydrological data using circulant embedding and dimensionality reduction. Water Resour. Res. 51, 4224–4243. https://doi.org/10. 1002/2014WR016395. Leeuwenburgh, O., Arts, R., 2014. Distance parameterization for efficient seismic history matching with the ensemble Kalman Filter. Comput. Geosci. 18 (3–4), 535–548. https://doi.org/10.1007/s10596-014-9434-y. (2014). Leonenkoy, Y., Keith, D., 2008. Reservoir engineering to accelerate the dissolution of CO2 stored in aquifers. Environ. Sci. Technol. 42, 2742–2747. Lindeberg, E., 1997. Escape of CO2 from aquifers. Energy Convers. Mgmt. 38, 235–240. Lindeberg, E., Zweigel, P., Bergmo, P., Ghaderi, A., Lothe, A., 2001. Prediction of CO2 distribution pattern improved by geology and reservoir simulation and verified by time lapse seismic. In Proceedings of the Fifth International Conference on Greenhouse Gas Control Technologies (GHGT-5) 299–304. Luo, X., Bhakta, T., Nædal, G., 2018. Correlation-based adaptive localization with applications to ensemble-based 4D seismic history matching. Spe J SPE-185936-PA. https://doi.org/l0.2ll8/l85936-PA. McLaughlin, D., Townley, L.R., 1996. A reassessment of the groundwater inverse problem. Water Resour. Res. 32 (5). https://doi.org/10.1029/2006WR005144. 1131e1161. Man, J., Zhang, J., Li, W., Zeng, L., Wu, L., 2016. Sequential ensemble-based optimal design for parameter estimation. Water Resour. Res. 52, 7577–7592. Metz, B., Davidson, O., de Coninck, H., Loos, M., Meyer, L., 2005. IPCC Special Report on Carbon Dioxide Capture and Storage Rep., Intergovernmental Panel on Climate Change. Working Group III, Geneva (Switzerland).

Acknowledgements This project was sponsored in part by support from Energi Simulation. The numerical simulations in this manuscript were performed using the Eclipse software donated by Schlumberger. The data and simulation results for this paper can be obtained by contacting the corresponding author at [email protected]. References Aanonsen, S.I., Nævdal, G., Oliver, D.S., Reynolds, A.C., Vallès, B., 2009. The ensemble Kalman filter in reservoir engineering–a review. Spe J. 14 (03), 393–412. Ajo-Franklin, J.B., Peterson, J., Doetsch, J., Daley, T.M., 2013. High-resolution characterization of a CO2 plume using crosswell seismic tomography: Cranfield, MS, USA. Int. J. Greenhouse Gas Control 18, 497–509. Alnes, H., Eiken, O., Nooner, S., Sasagawa, G., Stenvold, T., Zumberge, M., 2011. Results from Sleipner gravity monitoring: updated density and temperature distribution of the CO2 plume. 10th Int. Conf. Greenh. Gas Control Technol. 4, 5504–5511. Anchliya, A., Ehlig-Economides, C., Jafarpour, B., 2012. Aquifer management to accelerate CO2 dissolution and trapping, Paper SPE-126688-PA. Spe J. 17 (03). Bachu, S., Gunther, W.D., Perkins, E.H., 1994. Aquifer disposal of CO2: hydrodynamic and mineral trapping. Energy Convers. Manage. 35 (4), 269–279. Bauer, R.A., Jaques, P., Smith, V., Payne, W.G., 2014. General overview of IBDP microseismicity. Midwest Carbon Sequestration Science Conference 2014. Burnisona, S.A., Bossharta, N.W., Salakoa, O., Reed, S., Hamlinga, J.A., Goreckia, C.D., 2017. 4-D seismic monitoring of injected CO2 enhances geological interpretation, reservoir simulation, and production operations. Energy Procedia 114 (July), 2748–2759. https://doi.org/10.1016/j.egypro.2017.03.1539. Butsch, R., Brown, A.L., Bryans, B., Kolb, C., Hovorka, S., 2013. Integration of well-based subsurface monitoring technologies: lessons learned at SECARB study, Cranfield, MS. Int. J. Greenhouse Gas Control 18, 409–420. Cameron, D.A., Durlofsky, L.J., 2012. Optimization of well placement, CO2 injection rates, and brine cycling for geological carbon sequestration. Int. J. Greenhouse Gas Control. 10, 100–112. Cavanagh, A., 2013. Benchmark calibration and prediction of the Sleipner CO2 plume from 2006 to 2012. Energy Procedia 37, 3529–3545. Chadwick, R.A., Noy, D.J., 2010. History – matching flow simulations and time-lapse seismic data from the Sleipner CO2 plume. Vining, B.A., Pickering, S.C. (Eds.), Petroleum Geology: From Mature Basins to New Frontiers – Proceedings of the 7th Petroleum Geology Conference 1171–1182. Chen, Y., Zhang, D., 2006. Data assimilation for transient flow in geologic formations via ensemble Kalman filter. Adv. Water Resour. 29 (8), 1107–1122. Chen, Y., Oliver, D.S., 2017. Localization and regularization for iterative ensemble smoothers. Comput. Geosci. 21, 13–30. Cihan, A., Birkholzer, J.T., Bianchi, M., 2015. Optimal well placement and brine extraction for pressure management during CO2 sequestration. Int. J. Greenhouse Gas Control. 42, 175–187. Deutsch, C.V., Journel, A.G., 1998. GSLIB - Geostatistical Software Library and User’s Guide. Oxford University Press, New York - Oxford 1998. Doughty, C., Pruess, K., 2004. Modeling supercritical carbon dioxide injection in heterogeneous porous media. Vadose Zone J. 3, 837–847. ECLIPSE, 2012. ECLIPSE Reservoir Simulator (2012), Manual and Technical Description. Schlumberger GeoQuest, Houston. Emerick, A.A., Reynolds, A.C., 2013. History-Matching Production and Seismic Data in a Real Field Case Using the Ensemble Smoother With Multiple Data Assimilation. Society of Petroleum Engineershttps://doi.org/10.2118/163675-MS. Ennis-King, J.P., Paterson, L., 2003. Role of Convective Mixing in the Long-term Storage of Carbon Dioxide in Deep Saline Formations. Presented at Society of Petroleum Engineers Annual Technical Conference and Exhibition, Denver, Colorado SPE paper no. 84344. Espinet, A., Shoemaker, C., Doughty, C., 2013. Estimation of plume distribution for carbon sequestration using parameter estimation with limited monitoring data. Water Resour. Res. 49 (7), 4442–4464. Evensen, G., 1994. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res.: Oceans (1978–2012) 99 (C5), 10143–10162. Evensen, G., 2004. Sampling strategies and square root analysis schemes for the EnKF. Ocean Dyn. 54 (6), 539–560. Evensen, G., Hove, J., Meisingset, H., Reiso, E., Seim, K.S., Espelid, Ø., 2007. Using the EnKF for assisted history matching of a North Sea reservoir model. Paper Presented at SPE Reservoir Simulation Symposium.

214

International Journal of Greenhouse Gas Control 81 (2019) 199–215

W. Ma et al. Nævdal, G., Mannseth, T., Vefring, E.H., 2002. Near-well reservoir monitoring through ensemble Kalman filter. Paper presented at SPE/DOE Improved Oil Recovery Symposium. Oldenburg, C.M., Rinaldi, A.P., 2011. Buoyancy effects on upward brine displacement caused by CO2 injection. Transp. Porous Media 87 (2), 525–540. Oldenburg, C., Benson, S., 2002. CO2 injection for enhanced gas production and carbon sequestration. Society of Petroleum Engineers, SPE 74367, 2002. Oliver, D.S., Reynolds, A.C., Liu, N., 2008. Inverse Theory for Petroleum Reservoir Characterization and History Matching. Cambridge University Press. Pruess, K., Garcia, J., 2002. Multiphase flow dynamics during CO2 injection into saline aquifers. Environ. Geol. 42, 282–295 2002. Pruess, K., Bielinski, A., Ennis-King, J., Fabriol, R., Le Gallo, Y., Garcia, J., Jessen, K., Kovscek, T., Law, D.-S., Lichtner, P., Oldenburg, C., Pawar, R., Rutqvist, J., Steefel, C., Travis, B., Tsang, C.-F., White, S., Xu, T., 2003. Code intercomparison builds confidence in numerical models for geologic disposal of CO2. In: Gale, J., Kaya, Y. (Eds.), GHGT-6 Conference Proceedings: Greenhouse Gas Control Technologies, S, pp. 463–470. Roach, L.A.N., Angus, D.A., White, D.J., 2017. Assessment of the limitations on the seismic detectability of injected CO2 within a deep geological reservoir, 13th International Conference on Greenhouse Gas Control Technologies, GHGT-13, 14–18, November 2016, Lausanne, Switzerland. Energy Procedia 114, 4008–4019. https:// doi.org/10.1016/j.egypro.2017.03.1541. Sakurai, S., Hovorka, S.D., Ramakrishnan, T.S., Boyd, A., Mueller, N., 2005. Monitoring saturation changes for CO2 sequestration: petrophysical support of the Frio brine pilot experimen. SPWLA 46th Annual Logging Symposium. Society of Petrophysicists and Well Log Analysts. Sato, K., Mito, S., Horie, T., Ohkuma, H., Saito, H., Watanabe, J., Yoshimura, T., 2011. Monitoring and simulation studies for assessing macro- and meso-scale migration of CO2 sequestered in an onshore aquifer: experiences from the Nagaoka pilot site. Japan Int. J. Greenhouse Gas Control 5, 125–137. Shamshiri, H., Jafarpour, B., 2012. Controlled CO2 injection into heterogeneous geologic formations for improved solubility and residual trapping. Water Resour. Res doi: 10.1029 /2011WR010455, 48, W02530, 15 PP. Skjervheim, J., Evensen, G., Aanonsen, S.I., Ruud, B.O., Johansen, T.-A., 2007. Incorporating 4D seismic data in reservoir simulation models using ensemble Kalman filter. SPE J.-Richardson- 12 (3), 282. Spycher, N., Pruess, K., 2005. CO2-H2O mixtures in the geological sequestration of CO2.II. Partitioning in chloride brines at 12-100 C and up to 600 bar. Geochim. Cosmochim. Acta 69 (13), 3309–3320. Tarrahi, M., Jafarpour, B., 2012. Inference of permeability distribution from injectioninduced discrete microseismic events with kernel density estimation and ensemble

Kalman filter. Water Resour. Res. 48, W10506. https://doi.org/10.1029/ 2012WR011920. Tarrahi, M., Jafarpour, B., Ghassemi, A., 2015. Integration of microseismic monitoring data into coupled flow and geomechanical models with ensemble Kalman filter. Water Resour. Res. 51, 5177–5197. https://doi.org/10.1002/2014WR016264. Tian, L., Yang, Z., Fagerlund, F., Niemi, A., 2016. Effects of permeability heterogeneity on CO2 injectivity and storage efficiency coefficient. Greenhouse Gases Sci. Technol. 6 (1), 112–124. Tian, H., Pan, F., Xu, T., McPherson, B.J., Yue, G., Mandalaparty, P., 2014. Impacts of hydrological heterogeneities on caprock mineral alteration and containment of CO 2 in geological storage sites. Int. J. Greenhouse Gas Control. 24, 30–42. Trani, M., Wojnar, K., Moncorgé, A., Philippe, B., 2017. Ensemble-based assisted history matching using 4D seismic fluid front parameterization. Old Spe J. https://doi.org/ 10.2118/183901-MS. Van der Meer, L., 1992. Investigations regarding the storage of carbon dioxide in aquifers in the Netherlands. Energy Convers. Manage. 33 (5), 611–618. Van der Meer, L.G.H., Arts, R.J., Paterson, L., 2001. Prediction of migration of CO2 after injection in a saline aquifer: reservoir history matching of a 4D seismic image with a compositional gas/water model. In: Williams, D.J., Durie, B., McMullan, P., Paulson, C., Smith, A. (Eds.), Greenhouse Gas Control Technologies: Proceedings of the 5th International Conference on Greenhouse Gas Control Technologies, pp. 378–384. Vogt, C., Marquart, G., Kosack, C., Wolf, A., Clauser, C., 2012. Estimating the permeability distribution and its uncertainty at the EGS demonstration reservoir Soultzsous-Forêts using the ensemble Kalman filter. Water Resour. Res. 48https://doi.org/ 10.1029/2011WR011673. W08517. Vrugt, J.A., Robinson, B.A., 2007. Treatment of uncertainty using ensemble methods: comparison of sequential data assimilation and Bayesian model averaging. Water Resour. Res. 43 (1). https://doi.org/10.1029/2005WR004838. Wen, X., Chen, W., 2006. Real-time reservoir model updating using ensemble Kalman filter. Hesperia 11 (4), 431–442. White, D.J., 2013. Seismic characterization and time-lapse imaging during seven years of CO2 flood in the Weyburn field, Saskatchewan, Canada. Int. J. Greenhouse Gas Control. 16 (Suppl. 1), S78–S94. https://doi.org/10.1016/j.ijggc.2013.02.006. White, J.A., Chiaramonte, L., Ezzedine, S., Foxall, W., Hao, Y., Ramirez, A., McNab, W., 2014. Geomechanical behavior of the reservoir and caprock system at the In Salah CO2 storage project. Proc. Natl. Acad. Sci. U. S. A. 111, 8747–8752. Zhu, C., Zhang, G., Lu, P., Meng, L., Ji, X., 2015. Benchmark modeling of the Sleipner CO2 plume: calibration to seismic data for the uppermost layer and model sensitivity analysis. Int. J. Greenhouse Gas Control 43 (December), 233–246. https://doi.org/ 10.1016/j.ijggc.2014.12.016. 2015.

215