{'fraction_non_alphanumeric': 0.03289378350963629, 'fraction_numerical': 0.013911027387335169, 'mean_word_length': 5.010451140916217, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deep Learning-Accelerated 3D Carbon Storage Reservoir Pressure Forecasting Based on Data Assimilation Using Surface Displacement from InSAR

Fast forecasting of reservoir pressure distribution in geologic carbon storage (GCS) by assimilating monitoring data is a challenging problem. Due to high drilling cost, GCS projects usually have spatially sparse measurements from wells, leading to high uncertainties in reservoir pressure prediction. To address this challenge, we propose to use low-cost Interferometric Synthetic-Aperture Radar (InSAR) data as monitoring data to infer reservoir pressure build up. We develop a deep learning-accelerated workflow to assimilate surface displacement maps interpreted from InSAR and to forecast dynamic reservoir pressure. Employing an Ensemble Smoother Multiple Data Assimilation (ES-MDA) framework, the workflow updates three-dimensional (3D) geologic properties and predicts reservoir pressure with quantified uncertainties. We use a synthetic commercial-scale GCS model with bimodally distributed permeability and porosity to demonstrate the efficacy of the workflow. A two-step CNN-PCA approach is employed to parameterize the bimodal fields. The computational efficiency of the workflow is boosted by two residual U-Net based surrogate models for surface displacement and reservoir pressure predictions, respectively. The workflow can complete data assimilation and reservoir pressure forecasting in half an hour on a personal computer.


Introduction
Rapid data assimilation and associated reservoir performance forecasting is essential for robust management of many subsurface applications such as geologic carbon storage (GCS), geothermal energy recovery, oil/gas production and ground water resource management (Jiang & Durlofsky, 2020;Jo et al., 2021;Sun et al., 2021;Wu et al., 2021). The process of data assimilation (also known as stochastic inverse modeling) is to estimate the probability distributions of certain unknown parameters for a system by integrating observation data.
Forecasting the behavior of the system with quantified uncertainties is the chief goal of such data assimilation. Due to ubiquitous heterogeneity in subsurface porous media, we are usually concerned with a high-dimensional inverse problem that requires a large number of forward model runs to sufficiently capture the uncertainty. In the traditional workflow, these forward model runs and subsequent data assimilation are computationally expensive and time-intensive, which are major impediments to providing rapid data interpretations that can inform decisionmaking. Accelerating the data assimilation and the forecasting process will enable rapid decision-making regarding operational strategies and result in significant cost savings and more optimal system performance.
Estimating the pressure distribution in GCS reservoirs is important yet challenging. It is essential in the permitting phase as well as the monitoring, reporting, and verification (MRV) plan required under subpart of the US EPA Greenhouse Gas Reporting Program (Burton-Kelly et al., 2021;EPA, 2022). However, the spatially sparse observation data from a limited number of wells can lead to high uncertainties in the prediction of the reservoir pressure's spatial distribution. Accurate prediction of reservoir pressure distribution during CO2 injection requires high-dimensional observation data to track the reservoir pressure changes. Interferometric Synthetic-Aperture Radar (InSAR) has proven to be an effective monitoring technique to track the changes of reservoir surface deformation and thereby to infer reservoir pressure changes.
InSAR analysis compares a series of radar images of the same location, collected from orbiting satellites, at different times. Elevation changes of ground surface can be calculated from a detectable change in the phase of the reflected signal (Becker, 2006). InSAR has been widely used to measure ground surface displacements (Alghamdi et al., 2020;Hoffmann et al., 2001Hoffmann et al., , 2003. For GCS applications, InSAR has been successfully applied to measure surface uplift in the In Salah CO2 storage project, located in the central region of Algeria (White et al., 2014).
The data indicated an order of 5 mm per year ground surface uplift above the CO2 injection reservoir and an uplift pattern extending several kilometers around the injection wells (Vasco et al., 2008). The uplift magnitude can be explained by pressure-induced, poroelastic expansion of the storage reservoir based on coupled reservoir-geomechanics modeling (Rutqvist et al., 2010).
Those studies revealed that InSAR-measured deformation is highly correlated to reservoir pressure changes, a critical fact forming the basis of the current work. Compared with traditional well-based pressure monitoring, InSAR provides high-dimensional observation data at a low cost, potentially reducing the uncertainties in predicting the spatial distribution of reservoir pressure.
In this study, we proposed to develop a computationally efficient workflow to assimilate surface displacement data interpreted from InSAR to calibrate the underlying reservoir models and to predict dynamic reservoir pressure distributions. To obtain a desirable level of fidelity of reservoir pressure, the workflow is designed to calibrate heterogeneous rock properties such as rock facies, porosity, and permeability in three dimensions (3D). Most studies on data assimilation of surface displacement maps have been in two dimensions (2D) due to the high computational cost of 3D modeling. Iglesias & McLaughlin (2012) presented an inverse modeling approach to estimate heterogeneous permeability and elastic moduli in a 2D 40 × 40 gridded domain from measurements of surface displacement and production data. Recently, Alghamdi et al. (2020) developed a high-dimensional Bayesian inversion framework to infer permeability on a 16,896-cell 2D grid from surface deformation measurement. A full 3D data assimilation workflow requires a parameter space of at least 10 5 degrees of freedom. Such a high-dimensional parameter space renders traditional Bayesian approaches such as rejection sampling and Markov Chain Monte Carlo inapplicable. Ensemble-based methods such as ensemble smoother (ES) and ensemble Kalman filter (EnKF) are efficient alternative methods for high-dimensional inversion. Ensemble-based methods apply first-and second-order statistical information for error propagation, thus avoiding the high computational cost associated with gradient calculation (Ma et al., 2019;Tang, et al., 2021). In this study, we develop the workflow based on an iterative ES method called Ensemble Smoother with Multiple Data Assimilation (ES-MDA), developed by Emerick & Reynolds (2013). The ES-MDA approach has been demonstrated to have satisfactory performance with low computational cost Emerick & Reynolds, 2013;Liu et al., 2020).
To develop a computationally efficient workflow for rapid reservoir pressure forecasting, two major challenges must be addressed. The first challenge is to parameterize a highdimensional 3D parameter space which includes facies, porosity, and permeability distributions.
As explained in Section 4.1, we are concerned with bimodally distributed non-Gaussian fields represented by heterogeneous porosity/permeability distributions in each facies. Although deep neural networks such as generative adversarial networks (GAN) and variational autoencoders (VAE) have been recently applied in non-Gaussian field parameterization, studies that parameterize a 3D system with heterogeneous properties within individual facies are still scarce.
Two related works emerged recently and are both based on the successful application of convolutional neural networks. The first work applied a convolutional adversarial autoencoder network to generate an unconditional 3D non-Gaussian conductivity field (Mo et al., 2020). The second work is a variant of the principal component analysis (PCA) method, in which a convolutional neural network called "C3D net" is applied as a post-processor to PCA to improve the parameterization performance for a wide range of geological systems (Liu & Durlofsky, 2021). We adopted the parameterization method proposed in Liu & Durlofsky (2021), which has proven efficient and able to sufficiently honor the flow behavior of the system. The second challenge is to develop a computationally efficient surrogate model to predict surface displacement (a 2D map) from 3D non-Gaussian parameter distributions. The highdimensional and non-Gaussian nature of the problem makes traditional surrogate models such as Gaussian process inapplicable. Deep neural networks have been demonstrated as an effective solution for high-dimensional surrogate modeling by treating surrogate modeling as an image-toimage regression task. Mo et al. (2019) developed a deep dense convolutional neural network to predict CO2 plume migration in 2D Gaussian permeability fields. Identifying the large errors resulting when applying the model to 3D non-Gaussian conductivity fields, Mo et al. (2020) improved the model by adding multilevel residual blocks. Tang et al. (2021b) developed a 3D recurrent residual U-Net (R-U-Net) model to predict dynamic saturation and pressure fields for a binary channel system. The authors later applied the same neural network to predict dynamic saturation, pressure, and displacement fields in GCS applications with Gaussian permeability distributions (Tang et al., 2021a). In the current work, we develop a new R-U-Net architecture to predict 2D surface displacement maps directly from 3D non-Gaussian permeability and porosity distributions.
In summary, the major innovative contributions of this study include: 1) we develop a data assimilation workflow to predict high-resolution dynamic pressure distributions in GCS reservoirs using high-dimensional InSAR-based surface displacement data rather than sparse pressure measurements from wells; 2) we perform computationally efficient inversion for 3D non-Gaussian parameter fields; and 3) we develop a 3D-to-2D R-U-Net structure to directly predict surface displacement in the form of 2D maps.
The rest of the paper is organized as follows. First, we introduce the data assimilation and forecasting workflow with its major components. The forward numerical model and the relevant surrogate models applied in the workflow are introduced subsequently. We employ a commercial-scale 3D synthetic case to demonstrate the implementation of the workflow.
Numerical results for this case are presented and discussed. In the last section, we summarize our conclusions 2. Data Assimilation Workflow Fig.1 presents an overview of the data assimilation and forecasting workflow. There are four major components in the workflow: 1) model parameterization using the CNN-PCA approach, 2) data space dimension reduction using PCA, 3) ES-MDA inversion for model parameters based on InSAR data, and 4) surface displacement and reservoir pressure predictions with deeplearning surrogate models. We will describe the details of the first three components in the following sub-sections. The numerical forward model and the relevant deep learning surrogate models will be introduced in Section 3. Fig.1. Overview of the workflow for InSAR data assimilation and reservoir pressure forecasting. Each realization in the reservoir model ensemble is represented by a facies field, a porosity field, and a permeability field.

Parameterization of initial reservoir models
The purpose of model parameterization is to create a low-dimensional representation for initial reservoir models including facies, porosity, and logarithmic permeability distributions. To construct a traditional PCA representation, an ensemble of Nr realizations of vectorized geological variables ∈ ℝ ! ! ( " is the number of grid cells) is generated. These vectors are assembled into a data matrix ∈ ℝ ! " ×! # = 1 Where is the mean of the Nr realizations. A singular value decomposition yields = , where ∈ ℝ ! " ×! # and ∈ ℝ ! # ×! # are the left and right singular matrices, and ∈ ℝ ! # ×! # is a diagonal matrix containing singular values ( )) ) of . We select n largest components in and the corresponding columns in , denoted as * and * , respectively. A new PCA realization with features similar to that of original Nr realizations can be generated from any given lowdimensional vector l ∈ ℝ *×% sampled from (0, * ): l is also known as latent space vector. We can reconstruct the original realizations through the projected latent space vector (l ; ) ): The choice of the latent space vector dimension (n) can be informed by the criterion ≥ , where T represents the percentage of the total variance being preserved in the original data (Vo & Durlofsky, 2014). If n= Nr, = +,-,) will exactly match ) . This traditional PCA workflow is directly applicable to geological variables following a multi-Gaussian distribution. For non-Gaussian distributions (such as facies), the traditional PCA cannot be directly applied.
In this work, we will perform parameterization in two steps: the first step is to parameterize facies distributions using a CNN-PCA approach, and the second step is to parameterize porosity and permeability distributions using the traditional PCA approach. The CNN-PCA approach applies a transform network (fw) as a post processer for PCA models (Liu & Durlofsky, 2021): ,22.+,-is the resulting geological realization. For 3D models, a pre-trained C3D net (Tran et al., 2015) is applied as the transform network. The network is trained by minimizing the reconstruction error between +,-and = +,-,) . Additional training losses such as a style loss and a hard data loss are considered as well. The detailed formulations of each loss component can be found in the original paper.
At the second step, we construct a single traditional PCA model for porosity and logarithmic permeability distributions. This practice is based on a common situation where porosity and permeability are correlated, but the specific correlation between them is unknown.
For each realization of porosity and permeability distributions, ∈ ℝ '! ! is the concatenation of a normalized porosity vector and a normalized permeability vector, each with a length of " . The porosity and logarithmic permeability distributions are normalized based on their mean and standard deviation values in each facies of the Nr realizations.

Data parameterization
We use InSAR-inverted surface displacement map as our observation data. Each displacement map usually has a dimension of hundreds or thousands. Such a high-dimensional data space can lead to many problems such as extremely high computational cost and ensemble collapse issues in ensemble-based data assimilation methods. Therefore, we apply a similar traditional PCA approach as that presented in section 2.1 to reduce the dimension of the data space. To construct a PCA model, we form a data matrix ∈ ℝ ! & × (! ' 5%) combining the observation data vector 89: and predicted data vectors ; ( = 1, … , < ). The predicted data vectors are the surface displacement predictions based on the prior reservoir model realizations as shown in Fig.1. = denotes the data dimension, and Np denotes the number of the predicted data realizations.
= 1 = is the mean of the Np +1 data vectors. A singular value decomposition yields = . , and can be truncated accordingly to reduce the data space dimension from = to = ? ,. The latent space observation data I 89: ∈ ℝ ! ) * × % can be written as: Considering the forward model to be 89: = ( ) + , we can reformulate the problem in the latent data space as: where e is the measurement error based on the original data vector and = ? is the measurement error based on the projected data vector. Since ! ) * is a linear transformation matrix, the projected covariance matrix of the observation data measurement error ( ; = ) is correlated to the original error covariance matrix ( = ) through (Grana et al., 2019):

ES-MDA framework for model inversion
The ES-MDA approach updates model parameters (facies, porosity, and permeability) with observation data over several iterations to match the data to a satisfactory level. Here we introduce the implementation details of ES-MDA with an ensemble size of Ne.
• First, a prior ensemble of latent model parameters is randomly sampled from the standard normal distribution. We can generate the corresponding PCA models and prior reservoir model ensemble as described in section 2.1.
• Second, we use a surrogate model (to be introduced in section 3.2) to predict surface displacement ( ) ) based on the prior reservoir model ensemble, where = 1, 2, … , B . We can obtain the projected data vector ( I ) and I 89: ) and the projected covariance matrix of the data measurement error ; = ? following section 2.2.
• Third, the latent model parameter of the ith stochastic realization can be updated as follow: where = 1, 2, … F , denoting the iterations. l= ? C.% is the cross-covariance matrix between the model parameters and the projected data given by: l C.% and = ? C.% denote the mean of model parameter vector l ) C.% and data vector I ) C.% .
C.% is the auto-covariance matrix of the projected data given by: In Eq. (5), I DE,) C.% is the perturbed observation data (in the projected space) given by: . In our application, we set the of each iteration to % ! / .

Numerical Forward Model and Surrogate Models
The workflow requires generating responses, namely displacement and pressure distribution, from any given reservoir realizations. This can be done using either a high-fidelity numerical simulator or a neural network-based surrogate model. The latter is usually trained using data generated by the former. In the following sections, we will introduce the numerical forward model and the detailed neural network structure for the surrogate models.

High-fidelity forward model
The high-fidelity forward model couples multi-phase flow in porous media and geomechanics to predict reservoir pressure and surface displacement due to CO2 injection. The flow module describes the dynamics of CO2 being injected in a deep saline aquifer. We consider two components (CO2 and H2O) transported by two phases (the aqueous phase and the gaseous phase). The gaseous phase only contains the CO2 component and is in a supercritical state under subsurface conditions. The aqueous phase contains both CO2 and H2O components. The generalized mass conservation equation of component c is: where is porosity, denotes different phases, C is phase density, C is phase saturation, EC is the mass fraction of component c in phase l and C is the injected phase volume. C is phase velocity, which can be calculated based on Darcy's law: where k is the rock permeability, $C is the phase's relative permeability, C is phase viscosity, C is phase pressure, is the gravitational acceleration constant, and z is depth. The phase fractions and phase component fractions are calculated based on the model of Duan & Sun (2003). In the simulation, the phase densities and viscosities are updated via look-up in a precomputed table.
Surface displacement is calculated based on poromechanics (Biot, 1941). We solve the equilibrium equation of 3D porous media as follows: where ′ denotes effective stress, denotes the identity tensor, is the Biot coefficient, M is the density of porous media, and I is the equivalent pore pressure given by I = ∑ C C C . The effective stress tensor ′ is linked to the displacement field through a linear elastic constitutive relationship: is the rank-four stiffness tensor, is the solid displacement vector, and : is the inner product operator. In isotropic media, the stiffness tensor C (or Cijkl) can be calculated as: where )G is the Kronecker delta, K is the bulk modulus and G is the shear modulus. The pore pressure is first estimated by the flow module and then served as the input to the geomechanical module. The drained bulk modulus K is correlated to reservoir porosity through pore compressibility (cϕ): In this work, an open source multiphysics simulator GEOSX is used (Settgast et al., 2018).
Details of the numerical implementation and validation cases of the simulator can be found at http://www.geosx.org/. Although GEOSX can efficiently handle two-way full coupling, a oneway coupled approach is appropriate for our purpose: the pore pressure change and vertical effective stress change have a simple relationship as the vertical total stress remains largely constant.

Surrogate model based on 3D-to-2D R-U-Net
We develop two surrogate models in this section, one for surface displacement prediction and the other for depth-averaged reservoir pressure prediction. For both surrogate models, the output is a 2D map and the inputs are reservoir property distributions (3D fields) under a given injection scenario. We treat facies, porosity, and permeability fields in the storage reservoir as uncertain variables to be calibrated. Since the facies information is embedded in porosity and permeability distributions, we only use porosity and permeability distributions as inputs. Both surrogate models are based the same network structure as introduced below.

Application in a commercial-scale GCS model
The efficacy of the workflow is evaluated based on a synthetic commercial-scale GCS model.
The model has four injection wells and a two million metric tons total injection rate per year.
Model details and numerical results are introduced in the following sections.

Problem Setup
We build the reference reservoir model based on the clastic shelf depositional (AGD), which contains worldwide porosity/permeability measurements from deep saline aquifers. Fig.4(d) shows a cross plot between porosity and permeability based on each facies.
Like most of the field cases, there is not a simple, one-to-one correlation between porosity and permeability.
The four injection wells (represented by the thick black lines in Fig.4(a)) has a total injection rate of 2 million metric tons of CO2 per year. The four wells are grouped into a single manifold to honor the same tubing head pressure. Key reservoir properties are summarized in Table 2. The geomechanics model includes six overburden layers and one basement layer, of which the properties are summarized in Table 3. The reservoir is relatively thin compared with the vertical dimension of the formation. Note that all 3D visualizations of the reservoir in this paper are exaggerated by 100 times in the vertical dimension.

Geostatistical model and CNN-PCA realizations
To generate initial ensemble of reservoir models (i.e., informed guesses of reservoir property distributions), we use virtual measurements along nine well trajectories including four injectors (represented by the thick black lines in Fig.4(a)) and five exploration wells (represented by the thin black lines in Fig.4(a)). Virtual measurements are essentially queries of the reference reservoir model for rock facies, porosity, and permeability, which are referred as "hard data". We impose a 100% Gaussian noise in the measurement of porosity and permeability. Fig.5 shows the workflow to generate a single reservoir model realization from the hard data. Rock facies realizations are generated using Sequential Indicator Simulation (SIS). Normalized porosity and permeability realizations are generated using Sequential Gaussian Simulation (SGS) based on a selected random seed to preserve their dependency. Both SIS and SGS are conducted using the GSLIB software package (http://www.gslib.com/). The variogram parameters for SIS and SGS are summarized in Table 4. Both the variograms for SIS and SGS are assumed to be anisotropic with their major range and minor range sampled from a uniform distribution. After a rock facies realization is generated, porosity and permeability are rescaled based on the corresponding rock types.  We generate 1000 initial reservoir models from the workflow and present three realizations in Fig.6(a). We conduct PCA model parameterization as introduced in section 2.1 based on these 1000 initial reservoir models. The reduced dimension for facies ( % ) is 290 with 70% of total variance being preserved (T=70%). The reduced dimension for porosity and permeability ( ' ) is 656 with 90% of total variance being preserved (T = 90%). The total dimension of reduced parameter space is therefore 946. We reconstruct the reservoir model realizations in Fig.6(a) as shown in Fig.6(b). training process is to minimize the training loss (Ltrain) including a reconstruction loss and a weighted loss. The weighted loss is calculated at well locations. The mismatch between numerical model predictions ( ) ) and neural network predictions ( t ) is calculated using an L2 norm: where Y is the number of training reservoir model realizations, and Z is the number of time We apply the root mean square error (RMSE) and the coefficient of determination (R 2 ) as the metrics to evaluate model accuracy: where [ = 100 is the number of testing realizations. We test the accuracy of surrogate models trained with different numbers of training realizations, nm =1,000, 2,800, 4,500 and 6,200. Fig.7 presents the RMSE and R 2 metrics for displacement and pressure models. As nm increases from 1000 to 6,200, the displacement model accuracy increases from an R 2 value of 0.956 to 0.985, and pressure model accuracy increases from an R 2 value of 0.943 to 0.983. We apply the displacement model and the pressure model trained by 6,200 realizations in the following analyses and the workflow. respectively. Overall, the surrogate models have satisfactory predictive performance over a wide range of reservoir scenarios. The underperformance of realization 1 is mainly because the scenario where the upper right well takes almost no injection is not sufficiently represented in the training dataset. This scenario happens when the permeabilities along a wellbore are poorly connected vertically, which therefore reduces the injectivity of the well. We can also observe that the surface displacement distribution is highly correlated to the depth-averaged pressure distribution in each realization. Fig.8. (a) Surface displacement predicted by simulation (true) and the displacement surrogate model (prediction) after ten years of injection. (b) Depth-averaged reservoir pressure predicted by simulation (true) and the pressure surrogate model (prediction) after ten years of injection. Difference = prediction -true.

Data assimilation and forecasting results
We follow the workflow shown in Fig.1 to perform data assimilation and forecasting. We use an ensemble size B = 100 and an iteration number of F = 12 for ES-MDA. As mentioned in section 4.2, the parameter space dimension of this problem is 946 based on CNN-PCA parameterization. We consider surface vertical displacement at the ground surface after 2 years of injection as the observation data. We will use the observation data to calibrate the reservoir models to provide reservoir performance predictions at the observation time (2 years) and future times (after 2 years). A Gaussian white noise with a standard deviation of 5% of the observation values is applied to the simulated data to form the observation data as shown in Fig.10. We apply the data parameterization method introduced in section 2.2 to reduce the data space dimension.
In this application, > = 100 and = = 4096. With 90% of the total variance being preserved (T = 90%), the data space can be reduced to = ? = 11 latent variables for each realization.  To further examine whether the workflow can provide reasonable reservoir pressure predictions, we apply the pressure surrogate model to predict depth-averaged pressure maps after 2, 5 and 10 years of injection based on the posterior reservoir model assemble. As shown in Fig.12, the mean and standard deviation of each prediction ensemble is compared with the "true" pressure map from numerical simulation. We use RMSE to quantify the prediction performance.
In this application, t is the mean predicted pressure map from the posterior ensemble and ) is the "true" pressure map. The RMSE values are 0.08, 0.11, and 0.15 MPa respectively for three prediction periods. The results illustrate that the posterior reservoir models calibrated by surface displacement observation can provide reasonable predictions of reservoir pressure dynamics at current and future injection periods.   In this section, we analyze the effect of two factors on the performance of data assimilation and forecasting: observation data noise level and the degree of data space dimension reduction. As for the observation data noise level, we evaluate two Gaussian white noise levels: 7% and 10% of the observation values in addition to the baseline 5%. The other parameters of the workflow are kept the same as the previous analysis. We calculate the R 2 values of displacement and pressure predictions of each realization in the ensemble compared with the relevant true maps. Table 5 presents the 95% confidence interval of the calculated R 2 for injection durations of 2 and 10 years. The results show that increasing the Gaussian white noise level to 10% of the observation values will not significantly compromise the prediction performance. A higher noise level can lead to a wider uncertainty range of predictions. As for reduced data space dimension, we consider two additional preserved PCA variance percentages T = 94% and T = 97%, in addition to T = 90%. With Np = 100 being fixed, T = 94% results in a data space dimension of 15, and T = 97% results in a data space dimension of 24. A trial-and-error process reveals that we must increase the ES-MDA ensemble size (Ne) to avoid ensemble collapse issues as the data space becomes larger. We set Ne = 1000 for T = 94% and Ne = 3000 for T = 97%. For different preserved variance percentages, we summarize the 95% confidence interval of R 2 values for displacement and pressure predictions after 2 and 10 years of injection in Table 6. For both T = 94% and T = 97%, the prediction performance is slightly better than that of T = 90%. Preserving an appropriate percentage of variance is important: if the preserved variance percentage is too low, some key features of surface displacement maps might be missed, thereby undermining the fidelity of the inversion results.
On the other hand, if the preserved variance percentage is excessively high, some local features (such as the measurement noises) might be preserved to interfere with the inversion algorithm.  Table 7 summarizes the computational and training cost of the data assimilation and forecasting workflow. With an ensemble size of 100 and an iteration number of 12, the workflow can take approximately 3,000 core-hours to finish with GEOSX, a high-fidelity reservoir model, as the forward simulator. If we replace the forward simulator with the displacement and pressure surrogate models, the workflow runs in approximately half an hour on a 4-core laptop. The computational cost has been transferred to the generation of training datasets for the surrogate models, which costs approximately 14,000 core-hours. In a field application, the training process can be performed before data become available. The significantly accelerated data assimilation workflow allows multiple trials to determine the appropriate parameters for the workflow, which can significantly improve its efficiency in field applications. ~2 core-hours a 0.06 core-s b Single run (displacement) ~0.25 core-hours a 0.06 core-s b Data assimilation (100 ensemble, 12 iteration) ~2,925 core-hours a ~0.5 machine-hours b

Concluding Remarks
In this study, we develop a computationally efficient workflow to forecast reservoir pressure in 3D carbon storage reservoirs by assimilating surface displacement data. The surface displacement data is a 2D map that can be obtained from low-cost InSAR measurements. We employed an ensemble-based framework, ES-MDA, to drive the data assimilation. Various PCA techniques are applied to reduce the dimensionality of the parameter and data spaces. The significantly reduced parameter and data spaces boost the computational speed of the workflow by reducing the ensemble size required by ES-MDA. The workflow is also accelerated by two deep learning-based surrogate models for surface displacement and reservoir pressure predictions. The surrogate models are demonstrated to have good accuracy and generalizability to different durations of CO2 injection.
The workflow is evaluated against a synthetic commercial-scale GCS model with two facies. We apply a two-step PCA approach to tackle the parameterization challenge introduced by heterogeneous distributed permeability and porosity in individual facies. The same approach is also applicable to reservoir models with more than two facies. The workflow relies on additional well logs to provide hard data for the geostatistical model. We noticed for the synthetic case study, logs from the four injection wells alone are not sufficient to properly constrain the statistical models. Adding five virtual exploration wells successfully addressed this issue. This is not surprising as the concerned domain in the model is quite large, and a larger reservoir naturally requires more information to constrain. In practice, other data sources, such as a baseline 3D seismic inversion could be used to further constrain the models.
The efficacy of current workflow depends upon reliable surface displacement maps obtained from InSAR data inversion. A limitation of the current workflow is that the displacement maps we used were synthetic. Although a moderate amount of white noise was superposed, they perhaps are still "cleaner" than real InSAR data. The real-world application may involve additional challenges. For example, in a region with significant surface vegetation or other significant ground changes, the InSAR data between sweeps becomes decorrelated, leading to blocks of missing data. Correlated noise is also common due to some seasonal variations such as near-surface groundwater recharge. These are important "data engineering" issues to be addressed before the proposed workflow can be effectively applied to real-world applications. The contribution of the current work is mainly a "proof of concept"; a more robust solution accommodating the aforementioned issues would require more research and development effort in the future.