Fast uncertainty quantification of reservoir simulation with variational U-Net
Larry Jin Hannah Lu Gege Wen













Fast uncertainty quantification of reservoir simulation with variational U-Net
    Larry Jin,   Hannah Lu,  Gege Wen

Energy Resources Engineering

Stanford University

{zjin, hannahlu, gegewen}@stanford.edu
    December 30, 2023
================================================================================================================================





Quantification of uncertainty in production/injection forecasting is an important aspect of reservoir simulation studies. Conventional approaches include intrusive Galerkin-based methods (e.g., generalized polynomial chaos (gPC)  and stochastic collocation (SC) methods) and non-intrusive Monte Carlo (MC) based methods. Nevertheless, the quantification is conducted in reformulations of the underlying stochastic PDEs with fixed well controls. If one wants to take various well control plans into account, expensive computations need to be repeated for each well design independently. In this project, we take advantages of the equation-free spirit of convolutional neural network (CNN) to overcome this challenge and thus achieve the flexibility of efficient uncertainty quantification with various well controls. We are interested in the development of surrogate models for uncertainty quantification and propagation in reservoir simulations using a deep convolutional encoder-decoder network as an analogue to the  image-to-image regression tasks in computer science. First, a U-Net architecture is applied to replace conventional expensive deterministic PDE solver. Then we adopt the idea from shape-guided image generation using variational U-Net and design a new variational U-Net architecture for “control-guided" reservoir simulation. Backward propagation is learned in the network to extract the hidden physical quantities and then predict the future production by the learned forward propagation using the hidden variable with various well controls. Comparisons in computational efficiency are made between our proposed CNN approach and conventional MC approach. Significant improvements in computational speed with reasonable accuracy loss are observed in the numerical tests. 





§ INTRODUCTION



Reservoir simulation is an important tool for a wide range of engineering problems including nuclear waste management, oil and gas production, as well as carbon capture and storage (CCS). It helps us understanding these subsurface multiphase flow processes by solving spatially and temporally discretized mass and energy balance equations. However, the veracity of reservoir simulation results is uncertain due to the fact that the model parameters are subject to the uncertainties of various sources. Sources of parameter uncertainties include errors in data used for model parameterization, interpretive errors, and most importantly, heterogeneity of the subsurface under consideration. 

The goal of uncertainty quantification (UQ) in reservoir simulation is to evaluate the impact of parameter uncertainties on the model outputs. Conventional methods for UQ are mostly Monte Carlo (MC) based sampling <cit.>, which could induce huge computational costs in repetitive multi-scale/multi-physics deterministic simulations <cit.>. Moreover, well location, as another important control in reservoir simulation, needs to be considered in the evaluation of production/injection proposals. MC based method requires enough sampling in each fixed well location, which involves time-consuming computations. To mitigate the challenges in computation cost, and allow efficient quantification of uncertainty in flow dynamics with high dimensional random media, we introduce a CNN-based surrogate modeling tool (i.e., a variational U-net) in this project.


CNN-based generative adversarial network (GAN), autoregressive model, and auto-encoder (AE) are efficient tools to conduct image-to-image regressions. The recent applications of CNN-based models in reservoir simulation have achieved encouraging results <cit.>. In those models, the input reservoir properties (e.g., permeability) and output flow states (e.g., pressure and saturation) are treated as images.The CNN-based image-to-image regression models are therefore used as surrogate models for UQ tasks. However, the aforementioned frameworks are not able to incorporate information from controlling parameters such as bottom-hole pressure (BHP) or flow rates, which limit their potential in the UQ of reservoir models where control plays a distinct role. 

A variational U-net has been proposed by <cit.> in the context of image and video generation. A conditional U-net <cit.> for shape-guided image generation, conditioned on the output of a variational auto-encoder for appearance is presented, where the model learns to infer appearance from the queries and synthesizes images with specific appearances in different poses (i.e., shapes).

In this project, we propose to adopt both the U-Net and the variational U-Net as surrogate models for efficient UQ tasks in reservoir simulations. The shape-guided image generation process with variational U-Net is analogous to “control-guided" reservoir simulation, and can be intuitively adopted for our tasks. Some previous studies related to this project are summarized in section 2. In section 3, we set up the problem in a rigorous mathematical frame and explain the dataset in details. In section 4, we propose a new variational U-Net architecture for control-guided reservoir simulation. Motivations and detailed explorations are also illustrated in subsections. Numerical results and comparison with conventional MC methods are presented in section 5. Finally, we draw a conclusion of this project and discuss some future research directions in section 6.




§ RELATED WORK


The subsurface is such a complex system involving interactive physics in multiple scales that there is no accurate deterministic model. Thus probabilistic models have been explored to account for the uncertainties from model errors, model parametrization, heterogeneity of the environment and various geometries of the boundary, giving rise to extensive research interests in uncertainty quantification in reservoir engineering. A dominant strategy for such problems is to solve the deterministic problem at a finite large number of realizations of the random inputs using Monte Carlo sampling <cit.>. Variations of MC methods, including quasi-MC <cit.>, Multi-Level-MC(MLMC) <cit.>, stratified MC <cit.>, and etc, are designed to make more efficient sampling. Intrusive methods, (e.g., moment equations <cit.>, polynomial chaos <cit.>, stochastic collocation <cit.>, method of distributions <cit.>), as alternatives to MC simulations, have been widely studied in the past decades. Although efficient in some problems, intrusive methods, in general, require additional efforts in reformulation of the model and reconstruction of the deterministic solvers. In the meantime, it is well known that all intrusive methods suffer from curse of dimensionality (in random space). We refer to reviews on these topics in <cit.>.


Recently, deep learning has been explored as a competitive methodology across fields such as fluid mechanics <cit.>, hydrology <cit.>, bioinformatics <cit.>, high energy physics <cit.> and others. In particular, <cit.> adopted an end-to-end image-to image regression approach for surrogate modeling governed by stochastic PDEs with high-dimensional stochastic input in random porous media. In addition, the deep neural networks are set under a formal Bayesian formula  to enable the network to express its uncertainty on its predictions when using limited training data. <cit.> and our project both study two-dimensional, single phase, steady-state flow through a random permeability field. However, the emphasis of <cit.> lies in Bayesian deep learning dealing with high dimensional random inputs but with fixed single well control. In this work, we focus on uncertainty propagation with various well controls, which would allow existing well information to infer the incomplete/uncertain geological properties and evaluate potential well locations during simulation. In large-scale reservoir project like CCS, such quantification from simulations will be of significant importance and economic value in industry.


In CS231N class, we have learned two different approaches to image generation in the context of deep learning: Variational Auto-Encoder (VAE) <cit.> and Generative Adversarial Networks (GAN) <cit.>. In <cit.>, a conditional U-Net for shape-guided image generation, conditioned on the output of a VAE for appearance is presented. The separation between shape and appearance is carefully modelled and thus an explicit representation of the appearance, which can be combined with new poses, is obtained. Motivated by this work, we established the analogue between shape-guided image and control-guided saturation map. Similarly, the hidden appearance in image is analogous to the underlying permeability map. Under this framework, the variational U-Net is trained to learn the underlying Darcy's law. We take advantages of the equation-free nature of CNN, allowing the freedom of saturation map generation under various well controls. Efficient uncertainty quantification tasks can be conducted then to provide valuable evaluations of potential injection plans.




§ DATASET

In this project, we consider a carbon capture and storage (CCS) problem as the reservoir simulation model of interest, where the multiphase phase flow consists carbon dioxide (CO_2) as gas and water as fluid. CCS is an essential climate change mitigation technology to reduce the concentration of CO_2 in the atmosphere. CO_2 captured from concentrated sources, the atmosphere, or through bio-energy production, is compressed into a liquid and injected into deep geological formations for long term sequestration. <cit.>



 §.§ Problem statement

Let S be a image of CO_2 gas saturation map from a dataset 𝒮, which consists of different saturation samples from different injection well locations in various geological formation. The geological properties of interest are characterized by the permeability field 𝐤 and the injection well location is denoted by 𝐲. Physically, the relationship between the three quantities S,𝐲 and 𝐤 are characterized by the following Eq. <ref>:

    ϕ∂_t S+∇·𝐮 =0
       𝐮(𝐱) = -𝐤(𝐱,ω)∇Φ(𝐱), 𝐱∈𝒟,
       ∇·𝐮(𝐱) = f(𝐱),𝐱∈𝒟,
       𝐮(𝐱) ·𝐧̂(𝐱)=0 , 𝐱∈∂ D,
 
where ϕ is the porosity, 𝐮 is the velocity field and Φ is the potential including the pressure and gravity effect. The domain of interest (after projection) is 𝒟 = [0,l]×[0,w], where l is the length and w is the width of the system. The uncertainties of the system are indicated as the random variable ω in 𝐤(𝐱, ω). Without loss of generality, we consider a Gaussian permeability field 𝐤 with correlation length 25. In our simulation, the reservoir is rescaled in a 128× 128 mesh with 10-meter layer depth.  The spacial variable 𝐱 = (x_1,x_2) specifies the coordinates. f(𝐱) describes the injection by

    f(𝐱) = {   Q,     𝐱= 𝐲,
       0,    ,
        .

where Q is the injection rate.


We want to understand how the saturation maps S are influenced by the permeability map 𝐤 and well location map 𝐲 without repeatedly solving the equation system (<ref>). For an arbitrary given control of 𝐲, how can we predict the uncertainty propagation in the saturation map S will look like by using convolutional neural network? This can be interpreted as very useful information in the evaluation of potential CCS site.



 §.§ Data generation
 
The data set is generated by our team and the data preparation workflow is described as follows: 1) generate Gaussian permeability maps 𝐤; 2) prepare the random well location map 𝐲; 3) feed 𝐤 and 𝐲 to a numerical simulator and obtain the simulation results of CO_2 gas saturation map S and pressure map P; and 4) concatenate two sets of 𝐲, S, and P to make one data tuple.



In Figure <ref>, we demonstrate an example of the Gaussian permeability fields 𝐤 generated using the Stanford Geostatistical Modeling Software (SGeMS) <cit.>. The permeability varies from 0.001 mD to 200 mD to mimic the sandstone reservoir that can potentially be used for CCS projects.

The well location map 𝐲 incorporates the injection and the production information. There are a total of 10 injection wells on each well location map and each injection well has a constant supercritical CO_2 flow rate of 75 tons per day. Each map contains 6 production wells with a constant bottom hole pressure of 150 bar. An example of 𝐲 is shown in the first column of Figure <ref>. Note that the positive numbers of days indicate injection wells and the negative numbers of days indicate production wells.

The gas saturation maps S (second column in Figure <ref>) and pressure maps P (third column in Figure <ref>) are simulated with the state-of-the-art full-physics numerical simulator ECLIPSE <cit.>. Each simulation takes around 2 minutes on an Intel Core i7-4790 CPU.

Finally, we concatenate a sets of 𝐲, S, and P with a set of 𝐲', S', and P' (permuted) to become one data tuple as shown in Figure <ref>. The 𝐲, S, P, 𝐲', S', and P' within each tuple are always correlated to the same permeability map 𝐤 but different well location maps 𝐲. 




 §.§ Training, validation, and test sets
 
For the training and the validation set, a total of 42 different permeability map realizations 𝐤 (under the same Gaussian distribution) and 50 different well locations maps 𝐲 were simulated in 42 × 50 = 2100 runs. Each simulation lasts for 1,000 days and we took 4 time snapshots of the output, generating 2100 × 4 = 8400 sets of 𝐲, S, and P data and 8400 data tuples. Among the 8400 data tuples, we eliminated 16 tuples with extremely large pressure data, which are caused by numerical instability, and got a total of 8384 data tuples for training and validation. The training/validation set split is 7600/784 (roughly 10:1). 

For the test set, we considered a total of 100 different permeability map realizations 𝐤 (under the same Gaussian distribution as before) and 5 different well locations maps 𝐲, simulated in 100 × 5 = 500 runs. The tuple preparation process is similar to the training and the validation set. Since our test results is compared with a conventional MC method (more details in Section 4.1), the required number of permeability map scenarios are significantly higher than in the training set. In fact, rigorous conventional MC uncertainty quantification often requires more than 5,000 of different permeability map scenarios, which is computationally intractable in the scope of this project. The drawbacks of using a small number of permeability maps in the conventional MC is further discussed in Section 6.



§ METHODS

Our project is built up in two stages. In the first stage, we focus on replacing conventional PDE solver with neural network. For a deterministic system, the saturation map is predicted for a given pair of permeability field and well location using U-Net. Basically, it is supervised learning in the relationship of 𝐤,𝐲 and S governed by equation <ref>. Satisfactory results have be shown in section 5 of previous milestone report.

The second stage is to predict the saturation map for a given well location but without any prescribed knowledge of permeability. This is a more challenging and practical problem in CCS project. The detailed geological information of the field, e.g. the permeability, is always incomplete and uncertain. This project will be useful in the evaluation of the optimal well location even in the absence of the permeability information. In a real CCS project, data of the saturation map and pressure map will be collected from test well locations (usually not the optimal sites) in a pilot test. Geological information is then interpreted from the pilot test data to build up stochastic models. Then simulation is conducted to evaluate potential optimal well controls based on stochastic models. Our project aims to interpret the permeability from pilot test saturation and pressure map in existing test wells and predict the saturation maps in other potential well locations using the same deep learning framework. In other words, the network would be able to predict the saturation map in an arbitrary well location input different from the training sets. Efficient quantification can be conducted in the evaluation of optimal well controls.





 §.§ Conventional MC Review

We briefly review the conventional MC approach here and employ the MC solution as a baseline to compare with in the later experiment section. For each different well locations {𝐲_1,𝐲_2,⋯,𝐲_N}, one draw a set of permeability sample maps {𝐤_1(𝐱),𝐤_2(𝐱),⋯,𝐤_M(𝐱)} following the same Gaussian distribution. For these N× M pair (𝐲_i,𝐤_j), the following deterministic system is solved:

    ϕ∂_t S+∇·𝐮 =0
       𝐮(𝐱) = -k_j(𝐱)∇Φ(𝐱), 𝐱∈𝒟,
       ∇·𝐮(𝐱) = f(𝐱),𝐱∈𝒟,
       𝐮(𝐱) ·𝐧̂(𝐱)=0 , 𝐱∈∂ D,
       f(𝐱) = {   Q,     𝐱= 𝐲_i,
       0,    .
        .
 
The solution is recorded as S_i,j for different injection days: {S_i,j^t_1,S_i,j^t_2,⋯}. Notice that the saturation map under each well control is completely independent of each other. For a chosen well control 𝐲_i in a fixed injection day t_n, one can investigate quantities of interest from the sample set {S_i,1^t_n,S_i,2^t_n,⋯,S_i,M^t_n}. For example, the mean saturation map S̅_i^t_n can be approximated by

    S̅_i^t_n≈1/M∑_j=1^MS_i,j^t_n.

Law of large numbers guarantees the convergence of the above approximation with a rate of O(M^-1/2).





 §.§ Basic U-Net

The first step of our project is to do supervised learning using basic U-Net <cit.>. In <cit.>, a more appropriate network architecture, modified from the conventional U-Net <cit.>, is applied to single-phase flow problem. With the similarity of single-phase flow and our CO_2 injection, the U-Net architecture from <cit.> is adapted here. We simulate input images of {(𝐲_i,𝐤_i)} and the corresponding output target images {S_i} by equation (3). The architecture is illustrated in the following Figure <ref>:

Here z is the latent variable. The following regularized MSE loss is employed as the training loss function:

    L(𝐟_θ(𝐲,𝐤),S) = 1/n∑_i=1^n[𝐟_θ(𝐲_i, 𝐤_i)-S_i]^2+λθ_2,

where 𝐟_θ is the neural net depending on the parameter θ and λ is the regularization strength. We refer to Fig 1 in <cit.> for the detailed architecture of 𝐟_θ.

The performance of the prediction is quantified through the mean relative error, i.e.,

    e = 1/n∑_i=1^n𝐟_θ(𝐲_i,𝐤_i)-S_i_2^2/S_i_2^2,

where n is the total number of data points. Notice that everything is deterministic and explicitly forward (i.e., permeability information is needed as an input) at this point. The achievement of basic U-Net is to find a neural network that can replace the PDE solver for eqn (<ref>) efficiently. The results of Basic U-Net are presented in our milestone reports.



 §.§ Variational U-Net

Motivated by <cit.>, where a conditional U-Net for shape-guided image generation is combined with a VAE for appearance conditioned on the output, we find an analogy in <cit.> and our problem:

    ⟷
       ⟷
       ⟷


 Since the well location and the permeability capture all variations of saturation, the architecture proposed in <cit.> should be able to provide a saturation map generator controlled by the well location and permeability. We refer to Figure 2 and equation (3) in <cit.> for the architecture and loss function construction. However, the test results are not satisfactory. As we went through the loss function construction and the assumptions, we realized the potential cause of the problem. One key term in the ELBO of the loss function construction is the prior of appearance conditioned on shape, which captures potential interrelations between shape and appearance. However, there is no correlation between well location and the permeability field. These two variables are completely independent in our problem. Another issue is the assumption of Gaussian distribution of appearance conditioned on image and shape, which could be satisfied only if we change our previous training dataset.

Figure <ref> shows our idea of the modified architecture. For a pair of input data set (𝐲_i, S_i,j,P_i,j), which is saturation and pressure maps obtained from simulations in the existing old well locations 𝐲_i, autoencoder will learn the latent permeability 𝐤_j conditioned on  (𝐲_i, S_i,j,P_i,j). Here we involve the pressure map P_i,j, which can be conveniently obtained during the process, to better infer the underlying permeability. Then a U-Net will infer the new well location 𝐲_i' combined with the autoencoder so that the decoder in U-Net can generate a prediction of the new saturation map based on the new well location 𝐲_i' and the latent permeability 𝐤_j.



The loss function is modified as below:

    L(𝐟_θ,ϕ,ψ(𝐲_i'),S_i',j) 
            = 1/n∑_i=1^n[𝐟_ϕ(𝐲_i)-S_i,j]^2 
        +λ_θθ_2+λ_ϕϕ_2+λ_ψψ_2,

where λ_θ,λ_ϕ and λ_ψ are the regularization strength for θ, ϕ and ψ respectively. The performance of the predict, again is evaluated with the average relative error defined in Eq. <ref>. Notice that there the U-Net is combined with an AE not a VAE (as in <cit.>), and hence the loss function is defined as MSE loss here instead of the perceptual loss in <cit.>. The mean relative error of saturation (and similarly pressure) is defined the same as in eqn (6). The detailed implementation of E_θ, D_ψ and F_ϕ can be found in Appendix.

Compared with conventional MC approach, the benefits of the proposed method in computation efficiency come from two parts. First, the AE part substitute conventional PDE solver. Notice that the permeability information is interpreted from AE and not needed as an input anymore, which solves an inverse problem of importance already. Second, new well control is directly combined with underlying permeability in the network. The prediction will be saturation under a completely new well control that has never been solved by PDE or trained in the network. This enables freedom of well control without any additional computational costs (once training is complete). In the test process, one only needs to input saturation and pressure map from existing old well location, then the trained network will infer all the uncertainties of the underlying permeability and propagate the uncertainties to the saturation map in the new selected well locations. To achieve the same goal, traditional MC has no smart way but run expensive simulations again under different well locations in each permeability sample space. Hence, V-UNET enables us to get the uncertainty information in a fraction of time.

Quantity of statistic interests can be computed in the same fashion as conventional MC approach. For example, the mean saturation map S̅_i'^t_n in the new well location 𝐲_i' at injection day t_n can be approximated by

    S̅_i'^t_n≈1/M∑_j=1^M S_i',j^t_n.




  
 
   





  
   

 






§ EXPERIMENTS




The implementation of the variational U-Net following <cit.> is made with modifications on the model architecture. The model is trained with 7600 data points (train set). The validation set (dev. set) contains 784 data points. We test the results for UQ task on a separate data set of 500 data points (test set).


The quantities of interest in this context are the pressure and CO_2 saturation fields of the reservoir. As is mentioned in Section <ref>, permeability fields are only used for the full-physics numerical simulator to generate quantities of interest, and is not an input of our V-UNet model. We use mean relative error defined by Eq. <ref> as a metric of the performance. Note that the test error (shown in Table <ref> and <ref>) is the average of the mean relative error for both pressure and saturation.


We first run the experiments to predict the pressure and saturation fields for a single geological realization. A variety of model architectures and hyper-parameters are tested and compared. We also conducted the ablation study on the resNet layers between conv layers and at the low-dimension (see Appendix for detail) to see their impacts on the model. Finally, we compare our V-UNet based MC model with the conventional MC (baseline model) to find the overall speedup achieved and the accuracy loss. The results and discussion for each part are shown subsequently.


















 §.§ Single realization


We first run the experiments to predict the pressure and saturation fields for a single geological realization. A variety of model architectures and hyper-parameters are tested and compared, which are shown in Table <ref>. Note that the `normal' architecture in the table refers to the architecture demonstrated in Figure <ref>. The `reversed' architecture refers to the model that E_θ and F_ϕ in Figure <ref> are in reversed order (skip connections are established between F_ϕ and D_ψ). Models with varying learning rates (lr) and training epochs (ep) are also listed in the table.

We point out that one of the model settings does not give a converging loss value, so the associated test error is not shown. We can tell that the `normal' architecture with a learning rate of 2e-4 and epoch of 10 gives the lowest dev. loss, and the test error associated is 12.7%, which is a reasonably good accuracy for surrogate models in general. We will use this model setting as the `default' setting for the following ablation study.



A set of V-UNet results (with last set model architecture and hyper-parameters in Table <ref>) are shown in Figure <ref>. We show true and predicted saturation and pressure fields at two different time steps in the operation period. Row 1 and 2 show the saturation and pressure fields at 250 days, and row 3 and 4 show the saturation and pressure fields at 750 days. Column 1 for all rows shows the true solution from the data set, Column 2 shows the fields predicted by the V-UNet model, and Column 3 shows the difference map between Column 1 and 2.

As we can tell from the figures, the V-UNet predictions are in general in good agreement with the true fields, which clearly indicates the capability of V-UNet in capturing the saturation and pressure responses under a different set of well locations. The errors for saturation mostly comes from the edges of the CO_2 plumes, which could potentially be an indication of issues in predicting the underlying permeability fields or that the loss function is not sensitive to the location of the edge of the plume. The pressure field predictions are visually in close match with true solutions. However, we do see some `non-smoothness' in the predictions (more obvious in the difference maps), which indicates a violation of underlying PDEs. A physically-informed loss function (as is discussed in <cit.>) is a potential remedy of the issue.





 §.§ Ablation analysis


We perform an ablation analysis on the default V-UNet model, which comes from the last entry of Table <ref>. First, the resNet layers between each of the down-sampling conv layers in the encoder and between each of the up-sampling conv layers in the decoder are removed from the default model. Next, the resNet layers in the low dimension (the last three resNet layers in the encoder, and first three resNet layers in the decoder) are removed from the default model. The results of the ablation study is shown in Table <ref>. The ablated models give a higher dev. loss and test error compared to the default model, which indicates that both components are essential for the performance of the V-UNet. 






 §.§ UQ task compared to conventional MC


The V-UNet model is used as a surrogate model in the Monte-Carlo (MC) framework for the uncertainty quantification (UQ) task. The resulting V-UNet-based MC is compared to conventional MC (baseline model) to get a speedup versus accuracy loss. As a reminder, the goal of this task is that given the uncertainty information (e.g., simulation results with 100 different permeability realization) under a set of well locations, and we want to get the uncertainty for quantities of interest (e.g., pressure and saturation fields) under a different set of well locations. The conventional MC will have no other way but to re-run the simulations for new well locations, which is considerably time-consuming. However, the V-UNet will be able to infer the uncertainty results (with some accuracy loss) in a fraction of time.

The comparison of UQ results are shown in Table <ref>. As is mentioned earlier, the test data set contains 500 data points, which are 100 realizations for five of the well locations each. The average time in Table <ref> refers to the time required for simulation of all 100 realizations corresponding to one well locations and averaged over five well locations. The average error for pressure (P) and saturation (S) are calculated analogously.

Note that the baseline model (conventional MC) are used as reference solution, which means the errors on pressure and saturation fields are defined as 0.0%. The errors for V-UNet MC are calculated based on the difference to those of conventional MC. The error for pressure is 4.6% and that of the saturation is 21.84%, both of which are considered reasonable for surrogate models. The average time required for conventional MC is 12,200 seconds, while the that for V-UNet MC is 14 seconds, which corresponds to a 871 times speedup. We want to point out that a typical MC procedure will require at least 5000 simulation runs for a setting (we run 100 each due to the limitation on project time and computation resources), and running V-UNet model currently uses only a fraction of GPU memory. The speedup will be more significant in a realistic MC setting.



The mean saturation and pressure field for one of the five test well locations are shown in Figure <ref> and <ref>. In both sets of figures, the first sub-figure shows the true solution, the second sub-figure shows the predicted solution, while the third sub-figure shows the difference map of sub-figure 1 and 2. For both sets of figures, the fields predicted by V-UNet are in close agreement with the true solution, which demonstrate the possibility of using V-UNet as a surrogate model for fast uncertainty quantification of reservoir simulation with varying well locations.





We also want to point out that the V-UNet-based MC workflow significantly under predicts the variance of both pressure and saturation. Thus the results are not shown here. The issue is still to be investigated in the future work.





§ CONCLUSION AND FUTURE WORK


In this project, we explore the use of a variational U-NET to construct a surrogate model for a uncertain system governed by stochastic PDEs. The surrogate model is incorporated into the MC framework where quantities of statistic interest can be calculated. The image-to-image regression achieves satisfactory results in terms of predictive performance and uncertainty modeling. The main contribution of this work beyond previous studies is that we established a “control-guided" reservoir simulation, allowing variation of well control without any additional computational costs in an equation-free framework. In a MC-based uncertainty quantification framework, the proposed method accelerates computational speed significantly (compared with conventional MC) without much sacrifice of accuracy. CO_2 mean saturation can be predicted in new well locations with the overall plume feature captured, providing valuable evaluation for a large-scale CCS project.

There are several directions to improve and extend in our work. First of all, the MC sampling of 100 realizations is far away from statistically enough. Limited by the computational power, 100 sample of permeability does not provide enough uncertainty information in the geological space. This could explain our observations that 1. permeability is under represented compared to well location in our surrogate model; 2. the mean saturation is well-captured but the variance and higher statistic moments are not consistent with reference solution. In the future work, we need larger dataset to do more accurate uncertainty quantification (approximately 5000 permeability sampling). Alternatively, we can design a sampling in the low dimensional space incorporated in the network to give richer statistic variations as well as more speedup. A bayesian framework with corresponding perceptual loss needs to be considered as in the original V-NET model. In the regression part, physical-informed techniques <cit.> can be added to improve prediction accuracy on pressure. Last but not least, we plan to extend our work to time-dependent problem, where Recurrent Neural Network (RNN) could be incorporated in the current framework.













§ ACKNOWLEDGEMENTS



The experiments are conducted on the cluster of Stanford Center for Computational Earth & Environmental Sciences (CEES). The provided computational resources are greatly appreciated.



§ HONOR CODER CLARIFICATION

The Github source code corresponding to <cit.> (https://github.com/CompVis/vunet) is referred to during our implementation. We made significant modification with respect to the original code and wrote our version of V-UNet code from scratch (see source code in supplement material for detail).




§ APPENDIX



Here we show the detailed architecture of encoder and decoder models used in the V-UNet framework. The architecture of encoder (same layout but different model parameters for both E_θ and F_ϕ) is demonstrated in Table <ref>. Note that N_x = N_y = 128 denotes the number of x-y pixels in the input image. The variable n_v = 1 for E_θ since well location is the only type of input, while n_v = 3 for F_ϕ since we have saturation pressure and well location input at the same time.

The network-in-network layer (nin lyr) is simply a 1x1 convolutional layer with 8 filters. Both the `In res' and the `Low res' blocks denote the residual block with a structure of Conv2D-batchNorm-ReLU-Conv2D-batchNorm, where Conv2D stands for two-dimensional convolutional layer, batchNorm denotes batch normalization layer, and ReLU designates rectified linear unit. `In res' means the residual block in between each convolutional blocks, which is distinguished from `Low res' (i.e., the residual block with feature map in low dimensional space). The `Conv blk' refers to the convolutional block with a structure of conv2D-bathcNorm-ReLU.

Each of the outputs of the `In res' layer for E_θ will go through the skip connection of the U-Net and serve as the additional input (i.e., `skpInput') of the decoder model. The final outputs of both the encoder (with shape of 8×8×128 each), will be concatenated (shape of 8×8×256) and serve as the input of the decoder.



The detailed structure of the decoder (i.e, D_ψ) is shown in Table <ref>. The `low res' and `in res' are defined similarly as those in the encoder. The `upConv' stands for the upconvolutional block with a sequential stack of a two-dimensional unpooling layer, a reflection padding layer, a two-dimensional convolutional layer, a batch normalization layer, and a ReLU, which in the end expand the spatial dimension by two. The output of the upConv block is concatenated with the input from skip connection (`skpInput') with the exact same dimension before input to the next residual block. We have a conv2D layer in the end to convert the feature map to saturation and pressure fields.






ieee



