3D CNNPCA: A DeepLearningBased Parameterization for Complex Geomodels
Abstract
Geological parameterization enables the representation of geomodels in terms of a relatively small set of variables. Parameterization is therefore very useful in the context of data assimilation and uncertainty quantification. In this study, a deeplearningbased geological parameterization algorithm, CNNPCA, is developed for complex 3D geomodels. CNNPCA entails the use of convolutional neural networks as a postprocessor for the lowdimensional principal component analysis representation of a geomodel. The 3D treatments presented here differ somewhat from those used in the 2D CNNPCA procedure. Specifically, we introduce a new supervisedlearningbased reconstruction loss, which is used in combination with style loss and hard data loss. The style loss uses features extracted from a 3D CNN pretrained for video classification. The 3D CNNPCA algorithm is applied for the generation of conditional 3D realizations, defined on grids, for three geological scenarios (binary and bimodal channelized systems, and a threefacies channelleveemud system). CNNPCA realizations are shown to exhibit geological features that are visually consistent with reference models generated using objectbased methods. Statistics of flow responses (, , percentile results) for test sets of 3D CNNPCA models are shown to be in consistent agreement with those from reference geomodels. Lastly, CNNPCA is successfully applied for history matching with ESMDA for the bimodal channelized system.
Keywords:
Geological Parameterization Data Assimilation History Matching Deep Learning Principal Component Analysis Subsurface Flow∎
1 Introduction
In subsurface flow problems, data assimilation (history matching) typically entails the modification of geological parameters, such as permeability and porosity, such that flow simulation results are in essential agreement with observed data. Geological parameterization is often a useful component of this workflow, as it enables the mapping of highdimensional (gridblock level) descriptions to a set of uncorrelated lowdimensional variables. A key challenge for parameterization procedures is the preservation of the largescale geological structures that characterize the underlying geomodel. Recently, deeplearningbased parameterization algorithms have achieved promising results for the modeling of complex fluvial channel systems, though much of this work has been for 2D systems. Our goal in this study is to extend one such parameterization, CNNPCA (convolutional neural network – principal component analysis), to handle complex 3D geomodels.
Deeplearningbased treatments are now being applied for many aspects of subsurface characterization and flow modeling. This includes, e.g., surrogate modeling for flow (Mo2019; Wen2019; Tang2020; Jin2020), dataspace inversion for history matching and uncertainty quantification (Jiang2020), and dimension reduction of seismic data in timelapse seismic history matching (Liu2020b; Liu2020a). Initial work on deeplearningbased geological parameterization includes formulations involving fully connected autoencoders (AE) (Canchumuni2017) and deep belief networks (DBN) (Canchumuni2018). CNNbased algorithms are, in general, more efficient and scalable. Examples involving CNNs include convolutional variational autoencoders (VAE) (Laloy2017; Canchumuni2019a) and deep convolutional generative adversarial networks (GAN) (Chan2017; Chan2018; Dupont2018; Mosser2018; Laloy2018; Laloy2019; Chan2020; Azevedo2020). The combination of VAE and GAN has also been proposed by Mo2020 and Canchumun2020.
The abovementioned algorithms utilize deeplearning to establish an endtoend mapping between geological models and lowdimensional latent variables. The CNNPCA (Liu2019; Liu2020) and PCACycleGAN (Canchumun2020) algorithms take a different approach, in that deeplearning is used as a postprocessor for PCA, which is a traditional parameterization method. In 2D CNNPCA, a style loss based on geomodel features extracted from the VGG net (Simonyan2015a) is employed to improve geological realism. This idea has also been employed to improve the quality of parameterized models from a standard VAE (Canchumun2020). A wide range of deeplearningbased geological parameterization algorithms have been developed and tested on 2D systems. They have been shown to provide more realistic models than traditional methods such as standalone PCA (Sarma2006) and discrete cosine transform (Jafarpour2010). In a benchmark comparison performed by Canchumun2020, seven deeplearningbased parameterization techniques, including VAEbased, GANbased and PCAbased algorithms, were found to perform similarly for a 2D channel system.
The development and detailed testing of deeplearningbased parameterizations for 3D geomodels is much more limited. Most assessments involving 3D systems are limited to unconditional realizations of a single geological scenario (Laloy2018; Canchumuni2018; Mo2020). To our knowledge, conditional realizations of 3D geomodels have thus far only been considered by Laloy2017 for the Maules Creek Valley alluvial aquifer dataset (ti_lib). Thus there appears to be a need for work on deeplearningbased parameterizations for complex 3D geomodels.
In this study, we extend the 2D CNNPCA framework to 3D by incorporating several new treatments. These include the use of a supervisedlearningbased training loss and the replacement of the VGG net (appropriate for 2D models) with the C3D net (Tran2015), pretrained for video classification, as a 3D feature extractor. This formulation represents a significant extension of the preliminary treatment presented by TangLiu2020, where style loss was not considered (in that work geomodel parameterization was combined with surrogate models for flow in a binary system). Here we apply the extended 3D CNNPCA procedure for the parameterization of conditional realizations within three different geological settings. These include a binary fluvial channel system, a bimodal channel system, and a threefacies channelleveemud system. In all cases the underlying realizations, which provide the training sets for the 3D CNNPCA parameterizations, are generated using objectbased modeling within the Petrel geomodeling framework (manual2007petrel). Geomodels and results for flow statistics are presented for all three of the geological scenarios, and history matching results are presented for the bimodal channel case.
This paper proceeds as follows. In Section 2, we begin by briefly reviewing the existing (2D) CNNPCA algorithm. The 3D CNNPCA methodology is then discussed in detail. In Section 3, we apply 3D CNNPCA to parameterize facies models for a binary channel system and a three facies channelleveemud system. A twostep approach for treating bimodal systems is then presented. For all of these geological scenarios, 3D CNNPCA realizations, along with flow results for an ensemble of test cases, are presented. In Section 4, we present history matching for the bimodal channel system using the twostep parameterization and ESMDA. A summary and suggestions for future work are provided in Section 5. Details on the network architectures, along with additional geomodeling and flow results, are included in Supplementary Information (SI).
2 Convolutional Neural Networkbased Principal Component Analysis (CNNPCA)
In this section, we first give a brief overview of PCA and the 2D CNNPCA method. The 3D procedure is then introduced and described in detail.
2.1 PCA Representation
We let the vector , where is the number of cells or grid blocks, denote the set of geological variables (e.g., facies type in every cell) that characterize the geomodel. Parameterization techniques map to a new lowerdimensional variable , where is the reduced dimension. As discussed in detail in Liu2019, PCA applies linear mapping of onto a set of principal components. To construct a PCA representation, an ensemble of models is generated using a geomodeling tool such as Petrel (manual2007petrel). These models are assembled into a centered data matrix ,
(1) 
where represents realization , is the mean of the realizations, and the subscript ‘gm’ indicates that these realizations are generated using geomodeling software. A singular value decomposition of gives , where and are the left and right singular matrices and is a diagonal matrix containing singular values.
A new PCA model can be generated as follows,
(2) 
where and contain the leading left singular vectors and singular values, respectively. Ideally, it is the case that . By sampling each component of independently from the standard normal distribution and applying Eq. 2, we can generate new PCA models.
Besides generating new models, we can also apply PCA to approximately reconstruct realizations of the original models. We will see later that this is required for the supervisedlearningbased loss function used to train 3D CNNPCA. Specifically, we can project each of the realizations of onto the principal components via
(3) 
Here denotes lowdimensional variables obtained through projection. The ‘hat’ is added to differentiate from lowdimensional variables obtained through sampling (). We can then approximately reconstruct as
(4) 
where is referred to as a reconstructed PCA model. The larger the reduced dimension , the closer will be to . If all nonzero singular values are retained, will exactly match .
For systems where follows a multiGaussian distribution, the spatial correlation of can be fully characterized by twopoint correlations, i.e., covariance. For such systems, (constructed using Eq. 2) will essentially preserve the spatial correlations in , assuming is sufficiently large. For systems with complex geology, where follows a nonGaussian distribution, the spatial correlation of is characterized by multiplepoint statistics. In such cases, the spatial structure of can deviate significantly from that of , meaning the direct use of Eq. 2 is not appropriate.
2.2 2D CNNPCA Procedure
In CNNPCA, we postprocess PCA models using a deep convolutional neural network to achieve better correspondence with the underlying geomodels. This process can be represented as
(5) 
where denotes the model transform net, the subscript indicates the trainable parameters within the network, and is the resulting geomodel.
For 2D models, the training loss () for , for each training sample , includes a content loss () and a style loss ():
(6) 
where is a training set (of size ) of random new PCA models, and is a reference model (e.g., training image or an original realization ). Here quantifies the ‘closeness’ of to , and acts to ensure that the postprocessed model resembles (to some extent) the input PCA model. The style loss quantifies the resemblance of to in terms of spatial correlation structure.
As noted in Section 2.1, the spatial correlation of nonGaussian models can be characterized by highorder multiplepoint statistics. It is not practical, however, to compute such quantities directly. Therefore, in Eq. 6 is not based on highorder spatial statistics but rather on loworder statistics of features extracted from another pretrained CNN, referred to as the loss net (Gatys2015; Johnson2016). More specifically, we feed the 2D models , and through the loss net and extract intermediate feature matrices from different layers of the loss net. The uncentered covariance matrices, called Gram matrices, are given by , where and are the dimensions of . These matrices have been shown to provide an effective set of metrics for quantifying the multipoint correlation structure of 2D models. The style loss is thus based on the differences between and reference model in terms of their corresponding Gram matrices. This is expressed as
(7) 
The content loss is based on the difference between the feature matrices for and from a particular layer in the network.
For the 2D models considered in Liu2019 and Liu2020, VGG net was used as the loss net (Simonyan2015a). We have a concept of transfer learning here as the VGG net, pretrained on image classification, was shown to be effective at extracting features from 2D geomodels. However, since VGG net only accepts imagelike 2D input, it cannot be directly used for extracting features from 3D geological models. Our extension of CNNPCA to 3D involves two main components: the replacement of VGG net with 3D CNN models, and the use of a new loss term based on supervised learning. We now describe the 3D CNNPCA formulation.
2.3 3D CNNPCA Formulation
We experimented with several pretrained 3D CNNs for extracting features from 3D geomodels, which are required to compute style loss in 3D CNNPCA. These include VoxNet (Maturana2015) and LightNet (Zhi2017) for 3D object recognition, and C3D net (Tran2015) for classification of video clips in the sports1M dataset (Karpathy2014). These CNNs accept input as dense 3D tensors with either three spatial dimensions or two spatial dimensions and a temporal dimension. Thus all are compatible for use with 3D geomodels. After numerical experimentation with these various networks, we found the C3D net to perform the best, based on visual inspection of the geomodels generated by the trained model transform nets. Therefore, we use the C3D net in our 3D CNNPCA formulation.
We observed, however, that the Gram matrices extracted from the C3D net were not as effective as they were with VGG net in 2D. This is likely due to the higher dimensionality and larger degree of variability of the 3D geomodels relative to those in 2D. We therefore considered additional treatments, and found that the use of a new supervisedlearningbased loss term provides enhanced 3D CNNPCA geomodels. We now describe this procedure.
As discussed in Section 2.1, we can approximately reconstruct realizations of the original model with PCA using Eq. 4. There is however reconstruction error between the two sets of models. The supervised learning component entails training the model transform net to minimize an appropriately defined reconstruction error.
Recall that when the trained model transform net is used at test time (e.g., to generate new random models or to calibrate geomodels during history matching), this entails postprocessing new PCA models . Importantly, these new models involve rather than . In other words, at test time we do not have corresponding pairs of . Thus, during training, it is beneficial to partially ‘disrupt’ the direct correspondence that exists between each pair of . An effective way of accomplishing this is to perturb the reconstructed PCA models used in training. We proceed by adding random noise to the in Eq. 3; i.e.,
(8) 
where denotes the perturbed lowdimensional variable and is a perturbation term. Then we can approximately reconstruct with
(9) 
We now describe how we determine and specify . To find , we apply the ‘energy’ criterion described in Sarma2006 and Vo2014. This entails first determining the total energy , where are the singular values. The fraction of energy captured by the leading singular values is given by . Throughout this study, we determine such that the leading singular values explain ~ 80% of the total energy. For the components of the perturbation term , we set for , where is determined such that the first leading singular values explain ~ 40% of the total energy, and for . With this treatment, we perturb only the smallscale features in . This approach was found to be effective as it disrupts the precise correspondence between pairs, while maintaining the locations of major geological features.
We use the same realizations of as were used for constructing the PCA representation for the generation of . We reiterate that Eqs. 8 and 9 are used here (Eqs. 3 and 4 are not applied). The supervisedlearning loss function for each pair of , which we refer to as the reconstruction loss, is given by
(10) 
Note that in 3D CNNPCA we take in all cases.
The style loss is evaluated using a separate set of new PCA models , , with sampled from . As for the reference model , in 2D CNNPCA we used either a reference training image or one realization of . Here we generate realizations of the original geomodel using objectbased techniques in Petrel, so there is no reference training image. We therefore use realizations of to represent the reference. Instead of using one particular realization, all realizations of (in turn) are considered as reference models. Specifically, we use as the reference model for new PCA model . It is important to emphasize that and are completely unrelated in terms of the location of geological features – we are essentially assigning a random reference model () for each . However, because the style loss is based on summary spatial statistics, the exact location of geological features does not affect the evaluation of the loss.
The style loss between and the (noncorresponding) new PCA model is given by
(11) 
where are Gram matrices based on features extracted from different layers in the C3D net. The C3D net consists of four blocks of convolutional and pooling layers. Here we use the last convolutional layer of each block, which corresponds to . Details on the network architecture are provided in SI.
A hard data loss term is also include to assure hard data (e.g., facies type at well locations) are honored. Hard data loss is given by
(12) 
where is a selection vector, with indicating the presence of hard data at cell and the absence of hard data, and is the total number of hard data. The final training loss is a weighted combination of the reconstruction loss, style loss and hard data loss. For each pair of corresponding , and the unrelated new PCA model , the total training loss is thus
(13) 
The three weighting factors , and are determined heuristically by training the network with a range of values and selecting the combination that leads to the lowest mismatch in quantities of interest (here we consider flow statistics) relative to the original (Petrel) geomodels. We also require that at least 99.9% of the hard data are honored over the entire set of models. The training set is divided into multiple minibatches, and the total loss for each minibatch of samples is
(14) 
where is the batch size.
The model transform net for 3D CNNPCA is obtained by replacing the 2D convolutional layers, upsampling layers, downsampling layers and padding layers in 2D CNNPCA (Liu2019) with their 3D counterparts. This can be readily accomplished within the Pytorch deeplearning framework (Paszke2017). The training procedure for 3D CNNPCA is illustrated in Fig. 1. Each training sample consists of a pair of corresponding models and an unrelated new PCA model . The new PCA model and the reconstructed PCA model are fed through the model transform net . The reconstruction loss is evaluated using and the original model and Eq. 10. To evaluate the style loss, and are fed through the C3D net, and the relevant feature matrices are extracted. Then, Gram matrices are computed to form the style loss (Eq. 11). The final loss entails a weighted combination of reconstruction loss, style loss and hard data loss. The trainable parameters in are updated based on the gradient of the loss computed with backpropagation. This process iterates over all minibatches and continues for a specified number of epochs. The detailed architectures for the model transform and C3D nets are provided in SI.
3 Geomodels and Flow Results Using 3D CNNPCA
We now apply the 3D CNNPCA procedure to generate geomodels corresponding to three different geological scenarios (binary channelized, threefacies, and bimodal channelized systems). Visualizations of the CNNPCA models, along with results for key flow quantities, are presented. All flow simulations in this work are performed using Stanford’s Automatic Differentiation General Purpose Research Simulator, ADGPRS (zhou2012parallel).
3.1 Case 1 – Binary Channelized System
The first case involves a channelized system characterized by rock facies type, with 1 denoting highpermeability sandstone and 0 indicating lowpermeability mud. The geomodels are defined on a grid (144,000 total cells). The average sand fraction is 6.53%. Figure 2 displays four random facies realizations generated using objectbased modeling within Petrel (sandstone is shown in red, and mud in blue).



In this case, there are two production wells and two injection wells. The well locations are given in Table 1. All wells are assumed to be drilled through all 40 layers of the model, and hard data are specified (meaning ) in all blocks penetrated by a well. Wells are perforated (open to flow), however, only in blocks characterized by sand. These layers are indicated in Table 1.
A total of conditional realizations are generated to construct the PCA model (through application of Eq. 1). A total of singular values are retained, which explains 80% of the total energy. Then, reconstructed PCA models are generated using Eqs. 8 and 9. The first principal components explain 40% of the energy, so perturbation is only applied to for . A separate set of random PCA models is generated by sampling from .
P1  P2  I1  I2  

Areal location ()  (15, 57)  (45, 58)  (15, 2)  (45, 3) 
Perforated layers  18  24  1  8  15  22  1  8 
The realizations of , and form the training set for the training of the model transform net . The weighting factors for the training loss in Eq. 13 are , , and . These values were found to provide accurate flow statistics and nearperfect harddata honoring. The Adam optimizer (Kingma2014) is used for updating parameters in , with a default learning rate of and a batch size of . The model transform net is trained for 10 epochs, which requires around 0.5 hour on one Tesla V100 GPU.
After training, 200 new Petrel realizations and 200 new PCA models are generated. The new PCA models are fed through the trained model transform net to obtain the CNNPCA models. Truncation is performed on the new PCA models and on the CNNPCA models to render them strictly binary. Specifically, cutoff values are determined for each set of models such that the final sandfacies fractions match that of the Petrel models (6.53%). These models, along with the new Petrel realizations, comprise the test sets.
Figure 3 presents four testset PCA models (Fig. 3a), the corresponding truncatedPCA models (denoted TPCA, Fig. 3b) and the corresponding CNNPCA models (after truncation, Fig. 3c). Figure 4 displays the 3D channel geometry for a Petrel model (Fig. 4a), a truncatedPCA model (Fig. 4b), and the corresponding CNNPCA model (Fig. 4 c). From Figs. 3 and 4, it is apparent that the CNNPCA models preserve geological realism much better than the truncatedPCA models. More specifically, the CNNPCA models display intersecting channels of continuity, width, sinuosity and depth consistent with reference Petrel models.
In addition to visual inspection, it is important to assess the CNNPCA models quantitatively. We computed static channel connectivity metrics as well as flow responses for all of the test sets. The connectivity metrics suggested by PardoIguzquiza2003 were found to be informative in our 2D study (Liu2019), but for the 3D cases considered here they appear to be too global to capture key interconnectedchannel features that impact flow. Thus we focus on flow responses in the current assessment.
The flow set up involves aqueous and nonaqueous liquid phases. These can be viewed as NAPL and water in the context of an aquifer remediation project, or as oil and water in the context of oil production via water injection. Our terminology will correspond to the latter application.
Each grid block in the geomodel is of dimension 20 m in the and directions, and 5 m in the direction. Water viscosity is constant at 0.31 cp. Oil viscosity varies with pressure; it is 1.03 cp at a pressure of 325 bar. Relative permeability curves are shown in Fig. 5. The initial pressure of the reservoir (bottom layer) is 325 bar. The production and injection wells operate at constant bottomhole pressures (BHPs) of 300 bar and 340 bar, respectively. The simulation time frame is 1500 days. The permeability and porosity for grid blocks in channels (sand) are md and . Grid blocks in mud are characterized by md and .



Flow results are presented in terms of , and percentiles for each test set. These results are determined, at each time step, based on the results for all 200 models. Results at different times correspond, in general, to different geomodels within the test set.
Figure 6 displays fieldlevel results for the 200 testset Petrel models (red curves), truncatedPCA models (TPCA, blue curves), and CNNPCA models (black curves). Figure 7 presents individual well responses (the wells shown have the highest cumulative phase production or injection). The significant visual discrepancies between the truncatedPCA and Petrel geomodels, observed in Figs. 3 and 4, are reflected in the flow responses, where large deviations are evident. The CNNPCA models, by contrast, provide flow results in close agreement with the Petrel models, for both field and welllevel predictions. The CNNPCA models display slightly higher field water rates, which may be due to a minor overestimation of channel connectivity or width. The overall match in – results, however, demonstrates that the CNNPCA geomodels exhibit the same level of variability as the Petrel models. This feature is important for the proper quantification of uncertainty.
3.2 Impact of Style Loss on CNNPCA Geomodels
We now briefly illustrate the impact of style loss by comparing CNNPCA geomodels generated with and without this loss term. Comparisons of areal maps for different layers in different models are shown in Fig. 8. Maps for three PCA models appear in Fig. 8a, the corresponding TPCA models in Fig. 8b, CNNPCA models without style loss () in Fig. 8c, and CNNPCA models with style loss () in Fig. 8d. The flow statistics for CNNPCA models without style loss are presented in Fig. 9. Hard data loss is included in all cases. The CNNPCA models with reconstruction loss alone () are more realistic visually, and result in more accurate flow responses, than the truncatedPCA models. However, the inclusion of style loss clearly acts to improve the CNNPCA geomodels. This is evident in both the areal maps and the fieldlevel flow responses (compare Fig. 9 to Fig. 6 (right)).
We note finally that in a smaller 3D example considered in TangLiu2020, the use of reconstruction loss alone was sufficient to achieve welldefined channels in the CNNPCA geomodels (and accurate flow statistics). That case involved six wells (and thus more hard data) spaced closer together than in the current example. This suggests that style loss may be most important when hard data are limited, or do not act to strongly constrain channel geometry.
3.3 Case 2 – ThreeFacies ChannelLeveeMud System
This case involves three rock types, with the channel and levee facies displaying complex interconnected and overlapping geometries. The fluvial channel facies is of the highest permeability. The upper portion of each channel is surrounded by a levee facies, which is of intermediate permeability. The lowpermeability mud facies comprises the remainder of the system. The average volume fractions for the channel and levee facies are 8.70% and 5.42%, respectively. Figure 10 displays four realizations generated from Petrel (channel in red, levee in green, mud in blue). The width and thickness of the levees is 55–60% of the channel width and thickness. These models again contain cells. In this case there are three injection and three production wells. Well locations and hard data are summarized in Table 2.
We considered two different ways of encoding mud, levee and channel facies. The ‘natural’ approach is to encode mud as 0, levee as 1, and channel as 2. This treatment, however, has some disadvantages. Before truncation, CNNPCA models are not strictly discrete, and a transition zone exists between mud and channel. This transition region will be interpreted as levee after truncation, which in turn leads to levee surrounding channels on all sides. This facies arrangement deviates from the underlying geology, where levees only appear near the upper portion of channels. In addition, since we preserve levee volume fraction, the average levee width becomes significantly smaller than it should be. For these reasons we adopt an alternative encoding strategy. This entails representing mud as 0, levee as 2, and channel as 1. This leads to better preservation of the location and geometry of levees.
P1  P2  P3  I1  I2  I3  

Areal location ()  (48, 48)  (58, 31)  (32, 57)  (12, 12)  (28, 3)  (4, 26) 
Perforated layers in channel  15  21  25  30  1  5  15  20  25  30  1  6 
Perforated layers in levee      21  24    6  8  36  40 
We again generate Petrel models, reconstructed PCA models and a separate set of new PCA models for training. We retain leading singular values, which capture 80% of the total energy. The first 70 of these explain 40% of the energy, so perturbation is performed on , when generating the reconstructed PCA models. The training loss weighting factors in this case are , and . Other training parameters are the same as in Case 1.
Test sets of 200 new realizations are then generated. The testset PCA models are postprocessed with the model transform net to obtain the CNNPCA test set. The CNNPCA geomodels are then truncated to be strictly ternary, with cutoff values determined such that the average facies fractions match those from the Petrel models. Four testset CNNPCA realizations are shown in Fig. 11. These geomodels contain features consistent with those in the Petrel models (Fig. 10). In addition to channel geometry, the CNNPCA models also capture the geometry of the levees and their location relative to channels.
We now present flow results for this case. Permeability for channel, levee and mud are specified as 2000 md, 200 md and 20 md, respectively. Corresponding porosity values are 0.25, 0.15 and 0.05. Other simulation specifications are as in Case 1. Fieldwide and individual well flow statistics (for wells with the largest cumulative phase production/injection) are presented in Fig. 12. Consistent with the results for Case 1, we observe generally close agreement between flow predictions for CNNPCA and Petrel models. Again, the close matches in the – ranges indicate that the CNNPCA geomodels capture the inherent variability of the Petrel models.
TruncatedPCA geomodels, and comparisons between flow results for these models against Petrel models, are shown in Figs. S1 and S2 in SI. The truncatedPCA models appear less geologically realistic, and provide much less accurate flow predictions, than the CNNPCA geomodels.
3.4 Case 3 – Bimodal Channelized System
In the previous two cases, permeability and porosity within each facies were constant. We now consider a bimodal system, where logpermeability and porosity within the two facies follow Gaussian distributions. The facies model, well locations and perforations are the same as in Case 1. The logpermeability in the sand facies follows a Gaussian distribution with mean 6.7 and variance 0.2, while in the mud facies the mean is 3.5 and the variance is 0.19. The logpermeability at well locations in all layers is treated as hard data.
We use the sequential Gaussian simulation algorithm in Petrel to generate logpermeability values within each facies. The final logpermeability field is obtained using a cookiecutter approach. We use to denote the facies model, and and to denote logpermeability within the sand and mud facies. The logpermeability of grid block , in the bimodal system , is then
(15) 
Figure 13 shows four logpermeability realizations generated by Petrel. Both channel geometry and withinfacies heterogeneity are seen to vary between models.
For this bimodal system we apply a twostep approach. Specifically, CNNPCA is used for facies parameterization, and two separate PCA models are used to parameterize logpermeability within the two facies. Since the facies model is the same as in Case 1, we use the same CNNPCA model. Thus the reduced dimension for the facies model is . We then construct two PCA models to represent logpermeability in each facies. For these models we set , which explains ~ of the total energy. A smaller percentage is used here to limit the overall size of the lowdimensional variable (note that ), at the cost of discarding some amount of smallscale variation within each facies. This is expected to have a relatively minor impact on flow response.
To generate new PCA models, we sample , and separately from standard normal distributions. We then apply Eq. 2 to construct PCA models , and . The model transform net then maps to the CNNPCA facies model (i.e., ). After truncation, the cookiecutter approach is applied to provide the final CNNPCA bimodal logpermeability model:
(16) 
Figure 14 shows four of the resulting testset CNNPCA geomodels. Besides preserving channel geometry, the CNNPCA models also display largescale property variations within each facies. The CNNPCA models are smoother than the reference Petrel models because smallscale variations in logpermeability within each facies are not captured in the PCA representations, as noted earlier. The average histograms for the 200 testset Petrel and CNNPCA geomodels are shown in Fig. 15. The CNNPCA histogram is clearly bimodal, and in general correspondence with the Petrel histogram, though variance is underpredicted in the CNNPCA models. This is again because variations at the smallest scales have been neglected.
To construct flow models, we assign permeability and porosity, for block , as and . The median values for permeability within channel and mud facies are ~ 800 md and ~ 30 md, while those for porosity are ~ 0.16 and ~ 0.08. The simulation setup is otherwise the same as in Case 1.
Flow statistics for Case 3 are shown in Fig. 16. Consistent with the results for the previous cases, we observe close agreement between the , and predictions from CNNPCA and Petrel geomodels. There does appear to be a slight underestimation of variability in the CNNPCA models, however, which may result from the lack of smallscale property variation.
PCA geomodels and corresponding flow results for this case are shown in Figs. S3–S5 in SI. Again, these models lack the realism and flow accuracy evident in the CNNPCA geomodels.
4 History Matching using CNNPCA
CNNPCA is now applied for a history matching problem involving the bimodal channelized system. With the twostep CNNPCA approach, the bimodal logpermeability and porosity fields are represented with three lowdimensional variables: for the facies model, and and for logpermeability and porosity in sand and mud. Concatenating the three lowdimensional variables gives , with . This represents the uncertain variables considered during history matching.
Observed data include oil and water production rates at the two producers, and water injection rate at the two injectors, collected every 100 days for the first 500 days. This gives a total of observations. Standard deviations for error in observed data are 1%, with a minimum value of 2 m/day. The leftmost Petrel realization in Fig. 13 is used as the true geomodel. Observed data are generated by performing flow simulation with this model and then perturbing the simulated production and injection data consistent with standard deviations of measurement errors.
We use ESMDA (Emerick2013) for history matching. This algorithm has been used previously by Canchumuni2017; Canchumuni2018; Canchumuni2019a; Canchumun2020 for data assimilation with deeplearningbased geological parameterizations. ESMDA is an ensemblebased procedure that starts with an ensemble of prior uncertain variables. At each data assimilation step, uncertain variables are updated by assimilating simulated production data to observed data with inflated measurement errors. We use an ensemble size of . The prior ensemble consists of 200 random realizations of sampled from the standard normal distribution. Each realization of is then divided into , and . CNNPCA realizations of bimodal logpermeability and porosity are then generated and simulated.
In ESMDA the ensemble is updated through application of
(17) 
where represents simulated production data, denotes randomly perturbed observed data sampled from , where is a diagonal prior covariance matrix for the measurement error and is an error inflation factor. The matrix is the crosscovariance between and estimated by
(18) 
and is the autocovariance of estimated by
(19) 
In Eqs. 18 and 19, the overbar denotes the mean over the samples at the current iteration. The updated variables , , then represent a new ensemble. The process is applied multiple times, each time with a different inflation factor and a new random perturbation of the observed data. Here we assimilate data four times using inflation factors of 9.33, 7.0, 4.0 and 2.0, as suggested by Emerick2013.
Data assimilation results are shown in Fig. 17. The gray region denotes the – range for the prior models, while the dashed blue curves indicate the – posterior range. Red points show the observed data, and the red curves the true model response. We observe uncertainty reduction in all quantities over at least a portion of the time frame, with the observed and true data consistently falling within the – posterior range. Of particular interest is the fact that substantial uncertainty reduction is achieved in waterrate predictions even though none of the producers experiences water breakthrough in the history matching period.
Prior and posterior geomodels are shown in Fig. 18. In the true Petrel model (leftmost realization in Fig. 13), wells I1 and P1, and I2 and P2, are connected via channels. In the first two prior models (Fig. 18a, b), one or the other of these injector–producer connectivities does not exist, but it gets introduced in the corresponding posterior models (Fig. 18d, e). The third prior model (Fig. 18c) already displays the correct injector–producer connectivity, and this is indeed retained in the posterior model (Fig. 18f).
5 Concluding Remarks
In this work, the 3D CNNPCA algorithm, a deeplearningbased geological parameterization procedure, was developed to treat complex 3D geomodels. The method entails the use of a new supervisedlearningbased reconstruction loss and a new style loss based on features of 3D geomodels extracted from the C3D net, a 3D CNN pretrained for video classification. Hard data loss is also included. A twostep treatment for parameterizing bimodal (as opposed to binary) models, involving CNNPCA representation of facies combined with PCA for withinfacies variability, was also introduced.
The 3D CNNPCA algorithm was applied for the parameterization of realizations from three different geological scenarios. These include a binary fluvial channel system, a bimodal channelized system, and a threefacies channelleveemud system. Training required the construction of 3000 objectbased Petrel models, the corresponding PCA models, and a set of random PCA models. New PCA models were then fed through the trained model transform net to generate the 3D CNNPCA representation. The resulting geomodels were shown to exhibit geological features consistent with those in the reference models. Flow results for injection and production quantities, generated by simulating a test set of CNNPCA models, were found to be in close agreement with simulations using reference Petrel models. Enhanced accuracy relative to truncatedPCA models was also demonstrated. Finally, history matching was performed for the bimodal channel system, using the twostep approach and ESMDA. Significant uncertainty reduction was achieved, and the posterior models were shown to be geologically realistic.
There are a number of directions for future work in this general area. The models considered in this study were Cartesian and contained cells (144,000 total grid blocks). Practical subsurface flow models commonly contain more cells and may be defined on cornerpoint or unstructured grids. Extensions to parameterize systems of this type containing, e.g., cells, should be developed. Systems with larger numbers of wells should also be considered. It will additionally be of interest to extend our treatments to handle uncertainty in the geological scenario, with a single network trained to accomplish multiscenario style transfer. Finally, the development of parameterizations for discrete fracture systems should be addressed.
Computer Code Availability
Computer code and datasets will be available upon publication.
Acknowledgements.
We thank the industrial affiliates of the Stanford Smart Fields Consortium for financial support. We are also grateful to the Stanford Center for Computational Earth & Environmental Science for providing the computing resources used in this work. We also thank Obi Isebor, Wenyue Sun and Meng Tang for useful discussions, and Oleg Volkov for help with the ADGPRS software.References
Supplementary Information (SI)
In this SI we first provide details on the architectures of the 3D model transform net and C3D net. Modeling and flow results complementing those in the main text are then presented. These include truncatedPCA results for the threefacies example (Case 2) and PCA results for the bimodal channelized system (Case 3). Visualizations of geomodels and flow statistics are provided for both cases.
s.1 3D Model Transform Net Architecture
The architecture of the 3D model transform net is summarized in Table S1. In the table, ‘Conv’ denotes a 3D convolutional layer immediately followed by a 3D batch normalization layer and a ReLU nonlinear activation. ‘Residual block’ contains a stack of two convolutional layers, each with 128 filters of size and stride (1, 1, 1). Within each residual block, the first 3D convolutional layer is followed by a 3D batch normalization and a ReLU nonlinear activation. The second convolutional layer is followed only by a 3D batch normalization. The final output of the residual block is the sum of the input to the first 3D convolutional layer and the output from the second 3D convolutional layer. The last three ‘DeConv’ layers consist of a nearestneighbor upsampling layer, followed by a 3D convolutional layer with stride (1, 1, 1), a 3D batch normalization and a ReLU nonlinear activation. The last ‘DeConv’ layer is an exception as it only contains one 3D convolutional layer. Circular padding is used for all 3D convolutional layers.
Layer  Output size 

Input  (, , , ) 
Conv, 32 filters of size , stride (1, 1, 1)  (, , , ) 
Conv, 64 filters of size , stride (2, 2, 2)  (, , , ) 
Conv, 128 filters of size , stride (1, 2, 2)  (, , , ) 
Residual block, 128 filters  (, , , ) 
Residual block, 128 filters  (, , , ) 
Residual block, 128 filters  (, , , ) 
Residual block, 128 filters  (, , , ) 
Residual block, 128 filters  (, , , ) 
DeConv, 64 filters of size , upsample (1, , )  (, , , ) 
DeConv, 32 filters of size , upsample (, , )  (, , , ) 
Conv, 1 filter of size  (, , , ) 
s.2 C3D Net Architecture
The architecture of the C3D transform net is summarized in Table S2. In the C3D net, each ‘Conv layer’ contains a 3D convolutional layer followed immediately by a ReLU nonlinear activation. The C3D net contains four blocks of convolutional layers. The last convolutional layer at each block, corresponding to layers 1, 2, 4 and 6, is used for extracting features from 3D geomodels. These features are used for evaluating the style loss term. The 3D models are duplicated three times along the feature dimension (the fourth dimension) before being fed to the C3D net.
Layer  Output size 

Input  (, , , ) 
Conv layer 1, 64 filters of size  (, , , ) 
Maxpool layer, kernel size (1, 2, 2), stride (1, 2, 2)  (, , , ) 
Conv layer 2, 128 filters of size  (, , , ) 
Maxpool layer, kernel size (2, 2, 2), stride (2, 2, 2)  (, , , ) 
Conv layer 3, 256 filters of size  (, , , ) 
Conv layer 4, 256 filters of size  (, , , ) 
Maxpool layer, kernel size (2, 2, 2), stride (2, 2, 2)  (, , , ) 
Conv layer 5, 512 filters of size  (, , , ) 
Conv layer 6, 512 filters of size  (, , , ) 
Maxpool layer, kernel size (2, 2, 2), stride (2, 2, 2)  (, , , ) 
s.3 TruncatedPCA Results for ThreeFacies System
In this section we present truncatedPCA results for the threefacies (channelleveemud) system. For truncated PCA, the ‘natural’ encoding (mud denoted by 0, levee by 1, channel by 2) performs significantly better, both in terms of geomodel appearance and flow response, than the alternative encoding used for CNNPCA. Therefore, here we present truncatedPCA results using this natural encoding. See Section 3.3 in the main paper for a discussion of the alternate CNNPCA encoding.
The same set of Petrel realizations (encoded as noted above) is used to construct the PCA representation. The reduced dimension is , which captures ~ 80% of the total energy. These PCA realizations are truncated such that average facies fractions of mud, levee and channel match those of the Petrel models. Four testset truncatedPCA models are shown in Fig. S1. With reference to the Petrel realizations (Fig. 10 in the main paper), it is evident that the truncatedPCA geomodels do not display proper channel and levee geometries. Flow results for these models, presented in Fig. S2, exhibit significant deviation from the reference Petrel flow responses. The CNNPCA geomodels and flow results (shown in Figs. 11 and 12 in the main paper) are clearly superior to those shown here for truncated PCA.
s.4 PCA Results for Bimodal Channelized System
PCA results for Case 3 in the main paper will now be discussed. Here we construct the PCA representation directly for the bimodal logpermeability field . A total of Petrel realizations of the logpermeability field are used. For this case we set , which explains ~ 80% of the total energy. We also generate geomodels for which a histogram transformation is applied to the original PCA realizations, such that the average histogram for the resulting models matches that of the reference Petrel realizations. We refer to the PCA models after histogram transformation as PCA+HT models.
Two PCA models, along with the corresponding PCA+HT models, are shown in Fig. S3. By comparing these to the Petrel realizations (shown in Fig. 13 in the main paper), we see that PCA and PCA+HT models do not accurately reproduce the features evident in the reference geomodels. The PCA+HT models do, however, display better channel connectivity than the PCA models. The CNNPCA geomodels, shown in Fig. 14 in the main paper, more closely resemble the reference realizations.
Flow statistics for PCA and PCA+HT geomodels are presented in Figs. S4 and S5, respectively. Interestingly, the PCA models significantly underestimate flow quantities, while the PCA+HT models consistently overestimate these quantities. This may be due to a lack of channel connectivity in the PCA realizations, and overly wide channels in the PCA+HT models. Flow statistics for CNNPCA geomodels (Fig. 16 in the main paper), by contrast, display much closer overall agreement with the reference results from the Petrel models.