Integration of adversarial autoencoders with residual dense convolutional networks for estimation of non-Gaussian hydraulic conductivities
Shaoxing Mo Nicholas Zabaras Xiaoqing Shi Jichun Wu








Shaoxing Mo1,2, Nicholas Zabaras2, Xiaoqing Shi1, and Jichun Wu1

1Key Laboratory of Surficial Geochemistry of Ministry of Education, School of Earth Sciences and Engineering, Nanjing University, Nanjing, China.

2Center for Informatics and Computational Science, University of Notre Dame, Notre Dame, IN, USA.
 
Nicholas Zabarasnzabaras@gmail.com
Xiaoqing Shishixq@nju.edu.cn



    
  * A convolutional adversarial autoencoder is developed to parameterize non-Gaussian conductivity fields with multimodal distributions
    
  * A deep residual dense convolutional network is introduced as a surrogate of the forward physics-based model
    
  * The integrated method is tested with inverse problems for the estimation of non-Gaussian conductivities in solute transport modeling




Inverse modeling for the estimation of non-Gaussian hydraulic conductivity fields in subsurface flow and solute transport models remains a challenging problem.
This is mainly due to the non-Gaussian property, the non-linear physics, and the fact that many repeated evaluations of the forward model are often required. In this study, we develop a convolutional adversarial autoencoder (CAAE) to parameterize  non-Gaussian conductivity fields with heterogeneous conductivity within each facies using a low-dimensional latent representation. In addition, a deep residual dense convolutional network (DRDCN) is proposed for surrogate modeling of forward models with high-dimensional and highly-complex mappings. The two networks are both based on a multilevel residual learning architecture called residual-in-residual dense block. The multilevel residual learning strategy and the dense connection structure ease the training of deep networks, enabling us to efficiently build deeper networks that have an essentially increased capacity for approximating mappings of very high-complexity. The CCAE and DRDCN networks are incorporated into an iterative ensemble smoother to formulate an inversion framework. 
The numerical experiments performed using 2-D and 3-D solute transport models  illustrate the performance of the integrated method.
The obtained results indicate that the CAAE is a robust parameterization method for non-Gaussian conductivity fields with different heterogeneity patterns. The DRDCN is able to obtain accurate approximations of the forward models with high-dimensional and highly-complex mappings using relatively limited training data. The CAAE and DRDCN methods together significantly reduce the computation time required to achieve accurate inversion results.






§ INTRODUCTION


Groundwater flow and solute transport models are used widely to help understand subsurface processes and make science-informed decisions for groundwater resource management. 
Reliable model predictions that well reproduce the phenomena of interest require a good characterization of the hydraulic conductivity field as it greatly 
influences groundwater flow and solute transport. 
In many practical cases, such as aquifers in fluvial deposits where several highly contrasting facies coexist, it may be unrealistic to model the log-conductivity as a multivariate Gaussian. It has been shown that a multimodal 
distribution can better characterize the connectivity between different facies and the multimodality in highly heterogeneous conductivity fields <cit.>. 

In this study, we are concerned with an inverse estimation of non-Gaussian conductivity fields with heterogeneous conductivity within each facies using  hydraulic head and solute concentration data. This is in contrast to 
categorical fields with homogeneous conductivity within each facies (e.g., the binary field) <cit.>, although they both have multimodal distributions. In practice, the conductivity field is estimated via inverse modeling. Commonly used inverse methods are, for example, Markov chain Monte Carlo methods <cit.> and the ensemble-based data assimilation methods such as ensemble smoother <cit.>, ensemble Kalman filter <cit.>, and their variants <cit.>.
Considering the strong conductivity heterogeneity, the inverse problem 
is usually high-dimensional. Thus, it often requires a large number of forward model runs to obtain converged inversion results. 

To relieve the large computational cost, parameterization 
methods are commonly used together with surrogate models within the inversion framework. A parameterization method aims 
to represent the spatially correlated property field using a low-dimensional latent vector <cit.>. 
It can also 1) mitigate the potential ill-posedness of the inverse problem, 2) ensure that the updated fields in the 
inversion process satisfy the prior distribution assumptions (e.g., non-Gaussian) imposed on the unknown field, and 3) produce normally distributed latent variables, allowing direct use of second-order inverse methods (e.g., ensemble smoother and ensemble Kalman filter), although these may be at the cost of increased nonlinearity for the inverse problems <cit.>. In previous studies involving the estimation of non-Gaussian conductivity fields, as noted in the review by LINDE201586, one popular solution is to utilize multiple-point statistics (MPS) simulation in the inversion to generate conductivity field realizations that honor the  non-Gaussian prior determined by a training image <cit.>. While such MPS-based inverse methods can generally provide reliable estimations, they may be computationally intensive <cit.>. 
A surrogate method aims to replace the computationally expensive forward model with an accurate but cheap-to-run 
surrogate model during the inversion <cit.>. Although such combinations 
of methods for inverse modeling have been intensively studied 
for problems with Gaussian conductivity fields <cit.>, previous 
studies on problems with non-Gaussian conductivity fields often relied on the inverse methods solely without using the 
parameterization and surrogate methods. The development of parameterization and surrogate methods for such non-Gaussian problems remains an open problem due to two major challenges. 


First, most existing parameterization methods fail to work for non-Gaussian conductivity fields. Previous works on parameterizing the conductivity fields in inversion has relied on, for example, principal component analysis and its variants <cit.>. While these methods are well suited for Gaussian random fields, their performance for complex non-Gaussian fields deserves further improvement <cit.>. 
Inspired by the recent success of deep learning in various areas including Earth science <cit.> and hydrology <cit.>, its application in 
parameterization of non-Gaussian conductivity fields has been reported in many recent studies <cit.>. Among 
these applications, generative adversarial network (GAN) <cit.> and variational antoencoder (VAE) <cit.> 
are the two most popular network architectures. After training the network, these methods take random realizations of a low-dimensional vector as input and then 
generate new realizations of the conductivity field having similar features with those found in the training data. The quality of 
the generated realizations was shown to be superior to those from traditional parameterization methods <cit.>. 
However, these methods focused on  categorical conductivity fields with homogeneous conductivity within each facies. Their applicability to  non-Gaussian fields with heterogeneous conductivity within each facies, which is more challenging, remains to be further investigated. laloy2018 and Liu2019 further tested their methods' potential for such continuous fields. While the results suggested a promising performance of the deep learning-based parameterization methods, no inversion results were presented. Canchumuni201987 considered channelized conductivity fields with a bimodal distribution. They used a VAE network to generate binary channelized fields and the uncertain parameters to be estimated in the inversion are the latent variables and the permeability values at all grid points, with the inverse problem remaining high-dimensional.

Second, most existing surrogate methods suffer from the curse of dimensionality <cit.> and fail to efficiently obtain 
accurate approximations when the input-output relations are highly-nonlinear <cit.>. The curse of dimensionality is caused by the 
exponentially increased computational cost required for accurate surrogate construction as the input dimensionality (i.e., the number of uncertain variables considered) increases. 
Due to the strongly heterogeneous nature of the conductivity field, it is often required to use a large number of stochastic degrees of freedom to accurately represent the heterogeneity. 
The highly-nonlinear outputs here arise because the high-conductivity regions in a non-Gaussian aquifer result in preferential paths for the groundwater flow and solute transport.  The 
two factors together make the commonly used surrogate methods, such as Gaussian processes <cit.> and 
polynomial chaos expansion <cit.>, difficult to work. 
Deep neural networks have already exhibited a promising and impressive performance for surrogate modeling of forward 
models with high-dimensional input and output fields <cit.>. 
For example, in Tripathy2018 a deep neural network was proposed to build a surrogate model for a single-phase flow forward model. 
In Sun2018 and Zhong2019, their surrogate methods for a single-phase flow forward model and a multiphase flow forward model, respectively, 
were based on an adversarial network framework. In our previous studies <cit.>, 
a deep dense convolutional network (DDCN), which is based on a dense connection structure <cit.> for better information flow efficiency, was employed as the surrogate modeling framework. It showed a good performance 
in efficiently obtaining accurate surrogates of  various forward models 
with high-dimensional input-output mappings. However, these methods were tested on forward models with Gaussian 
conductivity fields <cit.> or on a single-phase 
flow model with binary channelized conductivity fields <cit.>. As shown in a case study below, the application of the DDCN surrogate method to solute transport modeling with non-Gaussian conductivity fields may lead to large approximation errors.  

In this work, we develop a convolutional adversarial autoencoder (CAAE) to parameterize  non-Gaussian conductivity fields with multimodal distributions. 
We transform a fully-connected adversarial autoencoder <cit.> to  a convolutional network so as 
to improve its scalability for larger-size inputs. 
In addition, we propose a deep residual dense convolutional network (DRDCN) for 
efficient surrogate modeling of forward models with high-dimensional and highly-complex mappings. Although deeper networks have the potential to 
substantially improve the network's performance, they can be difficult to train. We adopt in DRDCN a multilevel residual learning structure <cit.>. The residual learning strategy has been shown to be an effective solution to ease the training of very deep networks <cit.>. The multilevel residual learning structure is also implemented in the CAAE network. The CAAE and DRDCN networks are 
combined with an iterative local updating ensemble smoother (ILUES) algorithm <cit.> to formulate 
an efficient CAAE-DRDCN-ILUES inversion framework.  The overall integrated method 
is demonstrated using 2-D  and 3-D solute transport modeling with  non-Gaussian conductivity fields that have different heterogeneity patterns.  

In summary, three major innovative contributions are addressed in this study. First, we develop a CAAE method for parameterization of non-Gaussian conductivity fields with heterogeneous conductivity within each facies that is suitable in the context of inverse modeling. 
Second, we adopt a multilevel residual strategy in our previous DDCN method <cit.> to introduce a new DRDCN method with a substantially improved 
performance for surrogate modeling of highly-complex mappings. Finally and most importantly, to the best of our knowledge, we present the first attempt to incorporate simultaneously the parameterization and surrogate methods to perform inversion for non-Gaussian conductivities in solute transport modeling. 

The rest of the paper is organized as follows. In section <ref>, we introduce a solute transport model 
and define the problem of interest. The CAAE-DRDCN-ILUES inversion framework is presented in section <ref>. 
Then in sections <ref> and <ref>, the proposed method is evaluated using two synthetic examples. The conclusions are summarized in the last section.



§ PROBLEM DEFINITION

We consider  solute transport in heterogeneous porous media under a steady-state groundwater flow condition. It is 
assumed that the transport of solute is driven by  advection and dispersion. The governing equations 
for the steady-state flow and solute transport are written as <cit.>


    ∇·(K∇ h)=0,


and 


    ∂(ϕ c)/∂ t=∇·(ϕα∇ c)-∇·(ϕ vc)+r_s,


respectively. Here K (LT^-1) is the hydraulic conductivity, h (L) is the hydraulic head, ϕ (-) is the 
effective porosity, c (ML^-3) is the solute concentration, t (T) denotes time, r_s (ML^-3T^-1) 
is the sink/source, and α (L^2T^-1) is the dispersion 
tensor determined by the pore space flow velocity v (LT^-1), and longitudinal (α_L; L), transverse (α_T; L), and vertical (α_V; L) dispersivities. The two equations 
are coupled through the velocity v=-K/ϕ∇ h. The flow and solute transport equations are numerically solved using the MODFLOW <cit.> and MT3DMS <cit.> simulators, respectively. 

We are concerned with an inverse problem of characterizing the heterogeneous conductivity field using measurements of the hydraulic head and concentration.  The underlying conductivity fields of interest are non-Gaussian fields with a multimodal distribution. The inverse modeling is performed using the ILUES inversion algorithm <cit.> which 
has shown a promising performance for high-dimensional and highly-nonlinear inverse problems <cit.>.




§ METHODOLOGY




 §.§ ILUES for Inverse Modeling

The ILUES algorithm assimilates the output measurements d∈ℝ^N_d for multiple times with 
an inflated covariance matrix of the measurement errors to avoid overweighing the measurements <cit.>. 
The inflated covariance matrix is often taken as 𝐂̃_D=N_iter𝐂_D, where  𝐂_D  is 
the original covariance matrix of the measurement errors and N_iter is the number of iterations <cit.>. 
To better handle high-dimensional and highly-nonlinear problems, ILUES also adopts a local updating scheme, which updates each 
input sample m∈ℝ^N_m (in the CAAE-DRDCN-ILUES framework, m refers to the latent variables, see section <ref>) in the ensemble locally using its neighboring samples rather than all samples 
in the ensemble. Formally, given an ensemble of N_e input samples 𝐌=[ m_1,…, m_N_e], it first 
identifies a local ensemble for each sample m_i∈𝐌 based on the following metric <cit.>


    J( m)=J_d( m)/J_d^max + J_m( m)/J_m^max,


where J_d( m)=[f( m)- d]^⊤𝐂_D^-1[f( m)- d] quantifies the mismatch between 
the model responses f( m) and measurements d, and J_m( m)=( m- m_i)^⊤𝐂_MM^-1( m- m_i) 
is the distance between the sample m∈𝐌 and m_i. Here, 𝐂_MM is the autocovariance matrix 
of the input parameters in 𝐌, J_d^max and J_m^max are the maximum values of J_d(·) and J_m(·), 
respectively. Based on the J values, we select N_l=β_l N_e, (β_l∈(0,1]) samples as the local ensemble of m_i 
using a roulette wheel selection operator <cit.>, in which the selection probability 
of the ith individual is given as P_i=ρ_i/∑_j=1^N_eρ_j, i=1,…,N_e, where ρ_j=1/J( m_j) <cit.>. 
A local ensemble factor of β_l=0.1 suggested in zhang2018ILUES is used. 
 
 Let superscripts l, f, and a denote the local ensemble, current and updated samples, respectively. The ILUES first updates the local 
 ensemble of each sample m_i^f∈𝐌^f, that is, 𝐌_i^l,f, by using the usual ensemble smoother scheme <cit.>:

 
    m_i,j^a= m_i,j^f+𝐂_MD^l,f(𝐂_DD^l,f+𝐂̃_ D)^-1[d_j-f(m_i,j^f)],

 
 for j=1,…,N_l. Here 𝐂_MD^l,f is the cross-covariance matrix between 𝐌_i^l,f 
 and 𝐃_i^l,f=[f( m_i,1^f),…,f( m_i,N_l^f)], 𝐂_DD^l,f is the 
 autocovariance matrix of 𝐃_i^l,f, and d_j= d+𝐂̃_ D^1/2 r_N_d, r_N_d∼𝒩( 0,𝐈), 
 is the j-th realization of the measurements. The update of m_i^f, m_i^a, is then generated from its 
 updated local ensemble 𝐌_i^l,a=[ m_i,1^a,…, m_i,N_l^a] through a probabilistic scheme <cit.>. 
 One update iteration of ILUES is summarized in Algorithm <ref>. More details regarding ILUES 
 can be found in zhang2018ILUES and mo2019inverse.
 


 
For high-dimensional inverse problems, large ensemble size and iteration number are usually needed for ILUES 
to obtain converged and reliable inversion results, resulting in a large computational cost in forward model 
runs. To  reduce the computational burden, we propose a CAAE network for parameterizing 
the high-dimensional conductivity field using a low-dimensional latent vector and a DRDCN network 
to build an accurate but fast-to-run substitution of the forward model in the ILUES algorithm. 



 §.§ DRDCN for Surrogate Modeling

In the surrogate modeling task, we build a surrogate model to approximate the mapping between the input 
conductivity field and the output hydraulic head and concentration fields. In our previous 
studies <cit.>, we transformed the surrogate modeling task for problems 
with high-dimensional input and output fields  in a 2-D domain to an image-to-image regression problem by using a 
DDCN network. 
In this network, the input and output fields were treated as images. Denoting H× W as the 
spatial discretization resolution of the domain, 𝐱∈ℝ^n_x× H× W and 𝐲∈ℝ^n_y× H× W 
as the input and output fields, respectively. Then the surrogate modeling task for approximating the input-output mapping, 


    f: ℝ^n_x× H× W→ℝ^n_y× H× W,


was transformed to an image regression problem between n_x input images and n_y output images with a 
resolution of H× W, where n_x and n_y are the number of the input and output fields, respectively. 
It is straightforward to generalize to a 3-D domain by adding an extra depth axis to the images, that is, ℝ^n_x× D× H× W→ℝ^n_y× D× H× W. 

In order to further improve the performance of DDCN in problems with highly-complex mappings, we 
adopt a novel basic block called `residual-in-residual dense block' proposed in Wang-etal2018 for 
image super-resolution problems to formulate our DRDCN framework. 




  §.§.§ Residual-in-Residual Dense Block

A dense block introduces connections between non-adjacent layers aiming to fully exploit the hierarchical features 
from the outputs of preceding layers <cit.>. Let 𝐳^(i) (i=1,…,L) denote the output 
feature maps of the ith layer in the dense block, where L is number of layers. 𝐳^(i) is obtained by taking the concatenation of the output feature maps 
from its preceding layers as input, as represented by


    𝐳^(i)=ℋ([𝐳^(0),…,𝐳^(i-1)]),


where 𝐳^(0) represents the input to the dense block, and ℋ denotes operations on the input feature maps, including batch normalization (BN) <cit.>, 
followed by nonlinear activation and convolution (Conv) <cit.>. The ReLU (Rectified Linear Unit) function defined as (x)=max(0,x) is commonly used. Alternatively, a new activation function Mish <cit.>, which is defined as (x)=x·tanh(ln(1+e^x)), can be used. It has been shown to perform generally better than ReLU in many deep networks across many datasets <cit.>. 
A dense block with L=5 layers is illustrated in Figure <ref>a. 



It has been shown that deeper networks   have the potential to  
better approximate mappings of  high complexity, however,  they can be difficult to train <cit.>. 
To efficiently train a deeper network, we adopt   a multilevel residual learning structure, that is, the residual-in-residual 
dense block proposed in Wang-etal2018. In the residual learning framework, it has been shown that the residual mapping is much 
easier to learn than the original mapping <cit.>. Specifically, let f(𝐱) 
denote the desired underlying mapping to fit. Then the stacked layers (not necessarily the entire network) learn the residual mapping  
 g(𝐱):=f(𝐱)-𝐱.  The original mapping f(𝐱) is then 
recast as g(𝐱)+𝐱. Such a residual learning strategy can help alleviate the 
gradient vanishing problem for deep network training <cit.> and thus ease the training of very 
deep networks to achieve improved accuracy <cit.>. 

The architecture of residual-in-residual dense block is shown in Figure <ref>b. It consists of a stack of residual dense blocks, 
where the residual learning is used in two levels, resulting in a residual-in-residual structure. 
 That is, the residual learning implemented to the dense block results in a residual dense block; and that implemented to the stacked residual dense blocks results in a residual-in-residual dense block. 
The number of input and output feature maps of a dense block is N_ f (Figure <ref>a) and is set to N_ f=48. In addition to the residual-in-residual structure, a residual scaling technique <cit.> is also employed in the
residual-in-residual dense block to further increase the training stability <cit.>. Formally, this is done by scaling down the residual g(𝐱) 
by a factor β∈(0,1] before adding to 𝐱 (Figure <ref>b). A factor of β=0.2 suggested in Wang-etal2018 
is used in our network. 




  §.§.§ DRDCN Networks Based on Residual-in-Residual Dense Blocks

We employ the residual-in-residual dense block structure in our DRDCN network for surrogate modeling of solute transport in 
media with non-Gaussian conductivities. The network architecture is shown in Figure <ref>c. 
For 2-D or 3-D images (fields), the 2-D and 3-D operations (i.e., BN and Conv) implemented in the PyTorch software are respectively used in the network. The network contains 
four residual-in-residual dense blocks and the feature maps are to go through a coarsen-to-refine process. A convolutional layer 
is first employed to extract feature maps from the raw input image. The obtained features are then passed through the 
residual-in-residual dense blocks and the transition convolutional layers for downsampling/upsampling of the feature maps.
We arrange the position of the four residual-in-residual dense 
blocks in the network with a layout of (1,2,1) (Figure <ref>c) to encourage the information flow through the coarse feature maps. 
That is, one block is placed in the coarsening part; two adjacent blocks are placed in the most central part; and another one block is placed in the refining part.  
An additional level of residual learning is implemented to the stacked residual-in-residual dense blocks, resulting in a three-level residual learning structure in the network. In DRDCN, the Mish activation is used unless otherwise stated.





 §.§ CAAE for Parameterization of Non-Gaussian Random Fields

We propose to parameterize non-Gaussian conductivity fields with multimodal distributions using a CAAE network. Without loss of generality, here we use 𝐱 to denote the log-conductivity field. Adversarial 
autoencoder is a probabilistic autoencoder that uses the GAN framework as a variational inference algorithm <cit.>. 
The original adversarial autoencoder framework is composed of fully-connected layers <cit.>, making it increasingly 
difficult to train as the network gets deeper due to a large number of trainable parameters. To resolve this issue, we develop a 
CAAE framework based on convolutional layers to leverage their sparse-connectivity and parameter-sharing properties as well 
as robust capability in image-like data processing <cit.>. 




  §.§.§ Generative Adversarial Networks

GAN <cit.> is a framework that establishes an adversarial game between two networks: a generative network 𝒢(·) (i.e., generator) that learns the distribution p_data(𝐱) over the data, and a discriminative network 𝒟(·) (i.e., discriminator)  that computes the probability that a sample 𝐱 is sampled from 
p_data(𝐱), rather than generated by the generator. The generator maps the low-dimensional latent vector 𝐳 from 
the prior distribution p(𝐳) to the data space. The discriminator is trained to maximize the probability of 
distinguishing the real samples from the generated (fake) samples. The generator is simultaneously trained to maximally 
fool the discriminator into assigning a higher probability to the generated samples by leveraging the feedback from 
the discriminator. Mathematically, the adversarial game translates into the following minimization-maximization 
loss function <cit.>


    min _𝒢max _𝒟𝔼_𝐱∼ p_data(𝐱)[log𝒟(𝐱)]+𝔼_𝐳∼ p(𝐳){log[1-𝒟[𝒢(𝐳)]]}.


In practice, the generator and discriminator are usually trained in alternating steps: (1) train the discriminator to improve 
its discriminative capability; (2) train the generator to improve the quality of the generated samples so as to fool the discriminator.




  §.§.§ Adversarial Autoencoder

An autoencoder is a framework that learns a low-dimensional representation 𝐳 
(referred to as latent codes or latent variables) of a sample 𝐱 in the data and generates from the codes a 
reconstruction 𝐱̂ that closely matches 𝐱. It consists of two networks: an encoder to learn a 
mapping from 𝐱 to 𝐳 and a decoder to learn a mapping from 𝐳 to 𝐱̂. To create 
a generative framework, in adversarial autoencoder <cit.> a constraint is added on the encoder that forces it to generate latent codes 𝐳 
that roughly follow a desired distribution, like a standard normal distribution 𝒩(0,𝐈) 
used in the present study. The decoder is then trained to generate samples with features being consistent with those 
found in the training data given any sample 𝐳∼𝒩(0,𝐈) as input, resulting in a generative model. 

Mathematically, let q(𝐳|𝐱) be the encoding distribution, q(𝐱|𝐳) be  the 
decoding distribution, and p(𝐳) be the distribution that we want the latent variables 𝐳 to follow. The adversarial autoencoder looks for a generative model  


    p(𝐱)=∫ p(𝐱|𝐳)p(𝐳) d 𝐳.


The adversarial autoencoder <cit.> is very similar to the VAE <cit.> in the sense that 
in both a latent representation is obtained with a desired distribution. Thus we follow 
the formulation in VAE to finally introduce the adversarial autoencoder. In VAE, the generative model is obtained 
via minimizing the upper-bound of the negative log-likelihood


    𝔼_𝐱∼ p_data(𝐱)[-log p(𝐱)]< 𝔼_𝐱∼ p_data(𝐱){𝔼_𝐳∼ q(𝐳 | 𝐱)[-log p(𝐱 | 𝐳)]+KL[q(𝐳 | 𝐱)  p(𝐳)]
        }.


The first term on the right side quantifies the reconstruction quality and the second term is the Kullback-Leibler 
divergence measuring the difference between two distributions. 

In adversarial autoencoders, an adversarial training procedure instead of the Kullback-Leibler divergence is  used to
encourage an aggregated posterior distribution q(𝐳), instead of q(𝐳|𝐱), to 
match p(𝐳), where q(𝐳) is defined as <cit.>


    q(𝐳)=∫_𝐱 q(𝐳 | 𝐱) p_data(𝐱) d 𝐱.


An illustration of the adversarial autoencoder is depicted in Figure <ref>a. In the encoding path, the input 𝐱 is 
fed into the encoder which outputs two low-dimensional vectors of means μ and log-variances ln(σ^2) of the latent variables 𝐳. 
Then a vector 𝐳' is randomly drawn from 𝒩(0,𝐈) and rescaled to produce the codes 𝐳=μ+σ×𝐳', where × denotes element-wise 
multiplication.  The decoder takes 𝐳 as input to eventually generate 𝐱̂. Meanwhile, the discriminator of the 
adversarial network accepts input from the latent codes generated by the encoder or the prespecified 
distribution p(𝐳) to discriminatively predict whether the input arises from the encoding codes (fake sample) or p(𝐳) 
(real sample). Note that the adversarial network here differs slightly from the vanilla GAN framework <cit.>, 
in which the generator generates the sample 𝐱, and the discriminator discriminates 𝐱 whether it is from 
the generator or the data. 



The encoder (which is also the generator 𝒢(·) of the adversarial network), decoder, and discriminator 𝒟(·) of the 
adversarial autoencoder are trained jointly in two phases for each iteration: the reconstruction phase and the regularization phase <cit.>. 
In the reconstruction phase, the encoder (generator) and decoder are updated using the following loss function:


    ℒ_ED=ℒ_Rec+wℒ_𝒢,


where ℒ_Rec is the reconstruction error which in this study is taken as the L_1 loss:


    ℒ_Rec=1/N∑_i=1^N||𝐱_i-𝐱̂_̂î||_1,


and ℒ_𝒢 measures the generator's ability to fool the discriminator and has the form


    ℒ_𝒢=-1/N∑_i=1^Nlog{𝒟[𝒢(𝐱_i)]}.


Here, w is a weight factor balancing the two losses and a value of w=0.01 is used, 𝐱̂_i is the reconstruction of sample 𝐱_i, and N is the number of training samples. 
In the regularization phase, the discriminator is updated based on the loss function 


    ℒ_𝒟=-1/N∑_i=1^N{log[𝒟(𝐳_i)]+log[1-𝒟[𝒢(𝐱_i)]},


to distinguish the real sample from p(𝐳) (i.e. 𝐳_i) from the fake sample 𝒢(𝐱_i) produced by the generator. 

Such adversarial training process with the loss functions defined in equations (<ref>) and (<ref>) forces 𝐱̂ to 
closely match 𝐱 and q(𝐳) to gradually approach p(𝐳) (i.e., μ→0 
and σ^2𝐈→𝐈), respectively. After training, the decoder will define a generative 
model p(𝐱) that given an arbitrary input 𝐳∼𝒩(0,𝐈) can 
generate a new realization of sample 𝐱̂ with features similar to those in the data used for training (Figure <ref>d).




  §.§.§ CAAE Networks Based on Residual-in-Residual Dense Blocks

We also adopt the residual-in-residual dense block structure shown in Figure <ref>b in the encoder and decoder 
of the CAAE network. The encoder (Figure <ref>c) is similar to the coarsening part of the DRDCN network (Figure <ref>c). The encoder has two additional convolutional layers to respectively output the means μ and log-variances ln(σ^2). The decoder (Figure <ref>d) is similar to the refining part of the DRDCN network.  The decoder has an additional convolutional layer to 
extract feature maps from the codes 𝐳. The ReLU activation is used in the encoder and decoder. Inspired by Ledig2017CVPR, the discriminator is a stack of two convolutional layers followed by two fully-connected layers with 128 and 1 neurons, 
respectively. The leaky ReLU activation with a slope of 0.2 is used in the discriminator and the 
sigmoid activation is used in the last layer to output a probability value between 0 and 1. 




 §.§ The CAAE-DRDCN-ILUES inversion framework

We incorporate the CAAE parameterization method and the DRDCN surrogate method into ILUES to formulate an 
efficient inversion scheme for estimation of a non-Gaussian conductivity field of solute transport models. 
The integrated methodology is denoted as CAAE-DRDCN-ILUES hereinafter and is summarized in Algorithm <ref>. 
Notice that in this method, the surrogate model is used in Algorithms <ref> and <ref> 
to substitute the forward model. After parameterization, the uncertain parameters to be estimated are the latent variables 𝐳. 
The log-conductivity field is estimated with the following procedure: (1)  start with an initial latent code ensemble drawn 
from 𝒩(0,𝐈), (2) the corresponding log-conductivity fields are generated next using the CAAE's decoder, (3) the surrogate model is evaluated to obtain the predicted initial output ensemble, (4) the latent code and output ensembles are repeatedly updated for N_iter iterations using 
Algorithm <ref> based on the current latent code and output ensembles. In Algorithm <ref>, the input to the surrogate model to produce output predictions is the log-conductivity field, which is generated by the CAAE's decoder given the latent codes as input. The posterior log-conductivity fields are obtained from the decoder using the last latent code ensemble as input.







§ APPLICATION




 §.§ Solute Transport Models


The performance of the proposed method is illustrated using 2-D and 3-D solute transport modeling with random conductivity fields that have non-Gaussian heterogeneity patterns.



  §.§.§ 2-D Model

 The first test case considers 2-D solute transport within a channelized aquifer.  As shown in Figure <ref>a, the horizontal domain has a size of 10 (L)× 20 (L) and is uniformly discretized into H× W=40×80=3,200 gridblocks. The left and right 
boundaries are assumed to be constant head boundaries with heads of 1 (L) and 0 (L), respectively. 
No-flow boundary condition is imposed on the upper and lower boundaries. An instantaneous source with a concentration of 100 (ML^-3) is released from the location x=3 (L) and y=5 (L) at the initial time. The porosity and dispersivities 
are assumed to be known with constant values of ϕ=0.25, α_L=1.0 (L) and α_T=0.1 (L), respectively.



The non-Gaussian conductivity field in this case has a channelized and bimodal pattern with heterogeneous conductivity within each facies. The conductivity realizations are generated by the following procedure: First, 
a binary facies field is generated using the SNESIM code <cit.> with a training image shown in Figure <ref>c; then 
we populate each facies with log-conductivity values from two independently generated Gaussian random fields with a L_2 norm exponential covariance function:


    C(s,s')=σ_log(K)^2exp(-√((s_1-s'_1/λ_1)^2+(s_2-s'_2/λ_2)^2)),


where s=(s_1,s_2) and s'=(s'_1,s'_2) denote two arbitrary spatial locations, σ_log(K)^2 is the variance, and λ_1 
and λ_2 are the correlation lengths along the x- and y-axes, respectively. The means of Gaussian random fields corresponding to 
the high-conductivity channels and the low-conductivity non-channel medium are 4 and 0, respectively, while their variances and correlation 
lengths are taken the same with σ_log(K)^2=0.5, λ_1=4 (L), and λ_2=2 (L). 




  §.§.§ 3-D Model

The second test case considers  solute transport in a 3-D confined aquifer with a size of 1.5 (L)× 10 (L)× 20 (L) as depicted in Figure <ref>d. 
The domain is uniformly discretized into D× H× W=6× 32× 64=12,288 gridblocks. Similar to the 2-D case, the left and right boundaries are assumed to be 
constant head boundaries with heads of 1 (L) and 0 (L), respectively. No-flow boundary condition is imposed on the boundaries in the direction 
perpendicular to the y-axis. An instantaneous source with a concentration of 100 (ML^-3) is released from 
the location x=3 (L) and y=5 (L) at the initial time. The porosity and dispersivities are assumed to be 
known with constant values of ϕ=0.25, α_L=1.0 (L), α_T=0.1 (L), and α_V=0.01 (L), respectively. 
The conductivity fields for this 3-D model are obtained by randomly cropping 6× 32× 64 patches from a 120× 180× 150 training image depicted in Figure <ref>e (available at <http://www.trainingimages.org/training-images-library.html>). 
The conductivity heterogeneity pattern is different from that in the 2-D channelized field in the sense that the distribution of log-conductivities in this training image is trimodal with two major peaks around 1.0 and 3.0 and one minor peak around -0.6 (Figure <ref>g). 




 §.§ Synthetic Observations

As it will be shown in section <ref>, the CAAE-generated conductivity fields have higher regularity/smoothness 
than the original fields. In the inverse problem, the synthetic observations were obtained by running the forward model with 
an original conductivity field rather than a smoothed conductivity field generated by CAAE. The generation of synthetic observations 
here aims to mimic a real scenario, where data would be obtained by field measurements. The randomly generated 2-D and 3-D reference 
log-conductivity fields are depicted in Figures <ref>a and <ref>d, respectively. Note that they are distinct 
from the training log-conductivity samples of CAAE and DRDCN networks. The two reference fields are both non-Gaussian with a bimodal distribution (the 2-D case, Figure <ref>b) and a trimodal distribution (the 3-D case, Figure <ref>f), respectively. In the 2-D case, the 
concentration at t=[3,5,7,9,11] (T) and the hydraulic head are collected at 24 measurement locations (Figure <ref>a), resulting in 144 observations. 
In the 3-D case, the concentration at t=[4,6,8] (T) and the hydraulic head are collected at six depths of 24 measurement 
locations (Figure <ref>d), resulting in 576 observations. The synthetic observations were 
corrupted with 5% independent Gaussian random noise to the data generated by the reference model. Additionally, we do not use 
any conditioning data (i.e., measurements) of the conductivity, resulting in a rather challenging inverse problem. 




 §.§ Networks Design and Training




  §.§.§ CAAE Network Design and Training

The architecture of the CAAE network is shown in Figure <ref> and detailed in section <ref>. The encoder includes two downsampling layers which halves the feature map size via convolution with a stride of 2. Correspondingly, the two upsampling layers in the decoder double the feature 
map size to recover the output image size. The network consists of 71 layers, including 69 convolutional layers 
(mostly arising from the four residual-in-residual dense blocks that each contains 15 convolutional layers) 
and 2 fully-connected layers (in the discriminator). The kernel size in all convolutional layers is 3, and the stride in the convolutional 
layers that keep the same feature map size and halve the size is 1 and 2, respectively.

In the 2-D case, a training set with 40,000 realizations of the log-conductivity field is generated.

We also generate another 4,000 test realizations to evaluate the network's performance. In the 3-D case, we generate the log-conductivity realizations 
by cropping the 120× 180× 150 training image shown in Figure <ref>e. The original training image is split into two images, with the lower 105× 180× 150 part and the upper 15× 180× 150 part being used to generate the training and test datasets, respectively. The two images are both flipped 
along the three dimensions via the  operation implemented in MATLAB to augment the data, each resulting in four training images. The log-conductivity realizations with size 6× 32× 64 are then generated via cropping the training images using a stride of (2,6,10). In this way, we obtain 37,600 training samples and 3,000 test samples. The data augmentation strategy, which artificially creates new training data from existing training data via specific operations (e.g., flip, shift, and rotation), is commonly utilized to obtain improved performance <cit.>. Note that these operations may lead to correlation between the original and resulting samples, thus the training and test data are generated separately as mentioned above to avoid potential over-optimistic results when assessing the performance.
The reference log-conductivity fields used in section <ref> to generate the synthetic observations are randomly selected from the test sets.
The loss functions for network training are defined in equations (<ref>) and (<ref>). 
The network is trained on a NVIDIA GeForce GTX 1080 Ti X GPU for 50 epochs using the Adam optimizer <cit.> with a learning rate of 2× 10^-4 and a batch size of 64. The training required about 1.7 h and 13.1 h in the 2-D and 3-D cases, respectively.
 




  §.§.§ DRDCN Network Design and Training

The architecture of the DRDCN network is shown in Figure <ref> and detailed 
in section <ref>. The network is fully-convolutional and contains 64 convolutional layers without any fully-connected layers. 
The kernel size in all convolutional layers is 3, and the stride in the convolutional 
layers that keep the same feature map size and halve the size is 1 and 2, respectively. The softplus activation is used in the output layer for the concentration to ensure nonnegative predictions. Since the hydraulic head varies between 0 and 1 in both cases, the sigmoid activation is used in the output layer for the hydraulic head. 

The concentration at different time steps (i.e., t=[3,5,7,9,11] (T) and t=[4,6,8] (T) in the 2-D and 3-D cases, respectively) and the hydraulic head are collected as the observations in the inverse problem. 
Thus, the concentration fields at these time steps and the hydraulic head field are treated as the output channels of the network. 
There is one single input channel to the network which is the original log-conductivity field generated by following the 
procedure presented in section <ref>. In both cases, we generate four training sets with N=1,000, 2,000, 3,000, and 4,000 
samples to evaluate the convergence of the network approximation errors with respect to   the training sample size. 
The approximation accuracy is assessed using N_test=1,000 randomly generated test samples. Note that since the 3-D log-conductivity realizations are obtained via cropping the training image which may lead to data correlation issue (see section <ref>), the log-conductivity samples in the training and test data are randomly selected from the training and test datasets of the CAAE network, respectively. 
The accuracy is measured using the coefficient 
of determination (R^2), the root-mean-square error (RMSE), and the structural similarity index (SSIM) metrics. The R^2 and RMSE metrics are defined as
 
 
    R^2=1-∑_i=1^N_test||𝐲_i-𝐲̂_i||_2^2/∑_i=1^N_test||𝐲_i-𝐲̅||_2^2,

 
 and
 
 
    RMSE=√(1/N_test∑_i=1^N_test||𝐲_i-𝐲̂_i||_2^2),

 
 respectively, where 𝐲 denotes the simulated outputs, 𝐲̂ is the network predictions, and 𝐲̅=1/N_test∑_i=1^N_test𝐲_i. The SSIM metric measures the structural similarity between two images (fields) and is calculated over local windows of the image <cit.>


    SSIM(𝐮, 𝐯)=(2 μ_𝐮μ_𝐯+c_1)(2 σ_𝐮𝐯+c_2)/(μ_𝐮^2+μ_𝐯^2+c_1)(σ_𝐮^2+σ_𝐯^2+c_2),


where 𝐮 and 𝐯 are two windows with size 11×11 in the real and predicted images, respectively, μ_𝐮 (μ_𝐯) and σ_𝐮^2 (σ_𝐯^2) are the mean and variance values of window 𝐮 (𝐯), σ_𝐮𝐯 denotes the covariance between 𝐮 and 𝐯, and c_1=0.01 and c_2=0.03 are two constants <cit.>. Note that the SSIM metric is designed for 2-D images. A 3-D image with size D× H× W  is treated as D images with size H× W when computing the SSIM metric. A R^2 score value and a SSIM value approaching 1.0 and a lower RMSE value suggest better surrogate quality.
 
 
The network is trained using a regularized L_1 loss function:


    ℒ=1/N∑_i=1^N||𝐲_i-𝐲̂_̂î||_1+w_d/2θ^⊤θ,


where θ denotes all the network trainable parameters and w_d=1×10^-5 is a regularization coefficient. 
The network is trained on a NVIDIA GeForce GTX 1080 Ti X GPU for 200 epochs in the 2-D case and 300 epochs in the 3-D case using the Adam optimizer <cit.> with an initial learning rate of 
5× 10^-3 and a batch size of 32. We also use a learning rate scheduler which drops ten times on plateau during training. The training required about 0.3-1.2 h and 2.5-10.5 h in the 2-D and 3-D cases, 
respectively, as the training set size varies from 1,000 to 4,000.




§ RESULTS AND DISCUSSION

In this section, we first illustrate the performance of the CAAE and DRDCN networks in parameterization of non-Gaussian random fields 
and in surrogate modeling of the solute transport models, respectively. After that, the inversion results obtained 
from the CAAE-DRDCN-ILUES framework are compared to those obtained from the CAAE-ILUES framework without surrogate modeling.






 §.§ Parameterization of Non-Gaussian Random Fields

In the CAAE network, the latent dimensions in the 2-D and 3-D cases are set to 200 and 512, respectively. 
The parameterization results for  non-Gaussian log-conductivity fields 
are shown in Figure <ref>. It depicts the CAAE's reconstructions of three log-conductivity realizations in the test set, the corresponding histograms of the latent codes, 
and three random log-conductivity fields generated by the decoder with inputs from 𝒩(0,𝐈). 
In the reconstruction process, the original log-conductivity realization is fed to the encoder to produce latent codes, based on which the reconstructed field is generated by the decoder. 
It is observed that, in both 2-D and 3-D cases with different heterogeneity patterns, the network successfully recovers the spatial distributions of the low-conductivity and 
high-conductivity regions as well as the conductivities within these regions; although the conductivity heterogeneity is smoothed compared to the original fields. This heterogeneity smoothness is attributed to the information loss during the encoding and decoding processes.
The encoding latent codes roughly follow the prior distribution 𝒩(0,𝐈) 
that we imposed during training.  The reconstruction accuracy evaluated on the test sets is R^2=0.945 and 0.863 in the 2-D and 3-D cases, respectively. A larger latent dimension can be used in order to preserve more heterogeneity features in the generated fields, which however may lead to higher computational costs in inverse problems. The fourth row of Figures <ref>a and <ref>b 
shows the random log-conductivity fields generated by the decoder. The results show that the decoder is able to reproduce 
log-conductivity realizations that depict similar patterns of
heterogeneity (e.g., the channel structures and the conductivity continuity within the low/high-conductivity regions) to those found in the training data. 
Therefore, the CAAE network is employed in the inversion process as the parameterization framework for  non-Gaussian conductivity fields. 




 §.§ Surrogate Quality Assessment




To illustrate the superior performance of the proposed DRDCN network architecture against the DDCN network architecture 
employed in our previous studies <cit.> for surrogate modeling of systems 
with highly-complex input-output mappings, the DDCN network is also trained using the same training sets as those used in DRDCN. The DDCN network architecture is introduced in <ref>. 



The two networks' approximation accuracy for the 2-D and 3-D solute transport models is provided in Figure <ref>, 
which depicts the RMSEs, R^2 scores, and SSIMs evaluated on the test sets. It can be seen that DRDCN achieves lower RMSEs, higher R^2 
scores, and overall higher SSIMs than DDCN when being trained using the same four training sets in both cases. For example, with 4,000 training samples, our 
network achieves a RMSE of 0.016, a R^2 score of 0.972, and a median SSIM value of 0.994 in the 2-D case, while those obtained by the DDCN network are 0.024, 0.933, and 0.978, 
respectively (Figures <ref>a and <ref>c).  This implies that the DRDCN network can obtain accurate surrogates with fewer training samples 
(forward model runs) than the DDCN network. For instance, the 2-D DRDCN surrogate with 2,000 training samples results in  
approximation accuracy comparable to that of the DDCN surrogate with 5,000 training samples (thus with  50% reduction in training samples). The saved number of training samples indicates substantial 
computational gains especially for computationally intensive forward models in subsurface modeling where one single model execution 
can take up to hours or even days. The R^2 scores of the DRDCN networks for the concentration fields at different time steps and the hydraulic head field are depicted in Figures <ref>e and <ref>f. It is observed that the approximation accuracy for the hydraulic head is higher than that for the concentration. This is mainly because  the hydraulic head is less sensitive to the conductivity heterogeneity compared to the concentration <cit.>, leading to relatively smooth hydraulic fields that are easier-to-approximate.  In addition, there is no significant difference between the approximation accuracy for the concentration fields at different time steps. The DRDCN network's predictions for the output fields of 2-D and 3-D models given randomly 
selected input log-conductivity fields in the test sets are illustrated in Figures <ref> and <ref>, 
respectively. For comparison, the predictions by the forward models and network approximation errors are also shown in each plot. 
It is observed that although the output fields are highly irregular with sharp response changes, our network is able to obtain good approximations in both cases. 




It is worth noticing that the DRDCN network achieves higher accuracy improvement compared to the DDCN network 
in the 2-D case than in the 3-D case (Figure <ref>). The low-conductivity regions are barriers for groundwater flow and solute transport, leading to irregular output fields. In comparison to the 2-D log-conductivity fields which are mainly composed of low-conductivity regions (Figure <ref>c), the facies with the minimum average log-conductivity in the trimodal 3-D fields (i.e., the blue regions in Figure <ref>e) is only a small proportion of the aquifer as indicated in Figure <ref>g. The small proportion of low-conductivity regions in the 3-D fields thus results in relatively smoother output fields than those of the 2-D model. As a consequence, the DDCN network can also obtain a relatively accurate surrogate of the 3-D model, although its approximation error is larger than that of the DRDCN network. However, its performance greatly decreases in the 2-D model with more complex output fields. The results clearly suggest that the proposed DRDCN network performs better than the DDCN network in obtaining 
accurate surrogate results for the solute transport models with high-dimensional and highly-complex input-output relations. 
The improvement is attributed to the deeper network architecture with the help of the multilevel residual learning <cit.> 
and residual scaling strategies <cit.>. 
With a relatively small number of training samples, the DRDCN network is able to 
provide good approximations of the 2-D and 3-D forward models. Thus, it is used together with the CAAE parameterization 
strategy in the ILUES inverse method to formulate an efficient inversion framework for estimation of a non-Gaussian conductivity field of solute transport models. 




 §.§ Inversion Results


In both 2-D and 3-D cases, the DRDCN surrogates trained with 4,000 forward model runs are used to substitute the forward models in the CAAE-DRDCN-ILUES inversion framework. 
To assess the accuracy and computational efficiency of the surrogate-based method, the CAAE-ILUES method which evaluates the forward model rather than the surrogate model during the inversion is 
also performed.  The ensemble size of the ILUES algorithm is set to N_e=2,000 and N_e=3,000 in the 2-D and 3-D cases, respectively, to fully quantify the parametric uncertainty. 
The observations are assimilated for N_iter=30 iterations in the 2-D case and 20 iterations in the 3-D case. That is, the number of forward model runs required in CAAE-ILUES for the 2-D and 3-D cases are 31×2,000=62,000 (i.e., one prior ensemble and N_iter updated ensembles) 
and 21×3,000=63,000, respectively. The posterior log-conductivity fields are then obtained from 
the CAAE's decoder given the final latent code ensemble as input (Algorithm <ref>).  




The convergence of 
the fitting error between the model predictions and measurements as the iteration proceeds is shown 
in Figure <ref>, where the mismatch is measured using the normalized sum of squared weighted residuals (NSSWR): 


    NSSWR=1/SSWR_ref∑_i=1^N_d(f_i(𝐱̂)-d_i/σ_i)^2,


where {d_i}_i=1^N_d denote the measurements that contain measurement errors with standard deviations {σ_i}_i=1^N_d, and {f_i(𝐱̂)}_i=1^N_d are 
the forward model predictions given the input 𝐱̂ which is generated by the decoder. Here the SSWR metric is normalized using the reference value SSWR_ref. Thus, a NSSWR value approaching 1.0 suggests the convergence of the inversion process.
As discussed in sections <ref> and <ref>, the reference log-conductivity field has high heterogeneity versus the smoothed 
CAAE-parameterized field realizations 𝐱̂.  As a consequence, in the inversion process the NSSWR values are not able to converge to the reference NSSWR value (i.e., 1.0) due to the conductivity smoothness in the CAAE-generated samples. 
 This can be seen in Figure <ref> where 
although the NSSWR values in the ensemble approximately converge in both cases, 
the converged values (the mean values are 1.68 and 2.17 in the 2-D and 3-D cases, respectively) are larger than 1.0. Our tests showed that the NSSWR values were not able to converge asymptotically to 1.0 even when we increased the number of iterations in ILUES. The results of CAAE-DRDCN-ILUES are also shown in Figure <ref> (the surrogate prediction f̂(·) is used in equation (<ref>)). This method converges to larger NSSWR values (the mean values are 2.95 and 2.19 in the 2-D and 3-D cases, respectively) mainly due to the approximation errors of the DRDCN surrogate model. Accounting for the approximation errors in the inversion is expected to further improve the estimation accuracy <cit.>. This can be realized by building  a  model for the approximation errors and then accounting for the approximation error contribution to the covariance matrix 𝐂_D <cit.>. Alternatively, one can employ a Bayesian training strategy (e.g., Stein
variational gradient descent <cit.>) for the DRDCN network as in ZHU2018415, which can provide multiple sets of tuned network parameters and thus an estimation of the variance in the prediction. The uncertainty estimate can subsequently be incorporated in 𝐂_D as in zhang2016 to alleviate the influence of approximation errors. The training of Bayesian networks, however, is more expensive than that of non-Bayesian networks <cit.>.







Two posterior log-conductivity realizations, the ensemble mean, ensemble standard deviation, and histogram of the log-conductivities in the ensemble mean field of the CAAE-ILUES method for the two cases
are illustrated in Figure <ref>. The reference fields and the locations of the output 
measurements collected to infer the conductivity field are also shown in the plot to facilitate the comparison. 
It can be seen that in both cases the CAAE-ILUES can successfully capture the spatial distribution of the high-conductivity and low-conductivity regions as well as
the conductivities within these regions. This is also indicated by the multimodal distribution of the log-conductivities in the mean fields. Due to the relatively sparse measurements 
(i.e., only 24 observation wells are placed in the domain with thousands of gridblocks), the local conductivity 
heterogeneity, the location of the boundaries between high-conductivity and low-conductivity regions, the mean values of the modes in the histograms may not be accurately retrieved. 
As expected, the conductivity estimation in regions where no information is collected is less accurate with a larger estimation uncertainty 
than the estimation near the the output measurement locations. 
The results imply that the CAAE-ILUES algorithm performs well for this inverse problem but needs a 
large number of forward model evaluations. We show next the results 
of CAAE-DRDCN-ILUES in which the DRDCN surrogate models are used to substitute 
the 2-D and 3-D forward models. 
Various statistics for the same reference fields as in CAAE-ILUES 
are depicted in the third row of Figures <ref>a (2-D case) and <ref>b (3-D case).
Similarly, it can be seen that in both cases the surrogate-based 
framework successfully captured the multimodal features and identified the high-conductivity and low-conductivity regions as well as the conductivities within these regions. 
 





Figures <ref> and <ref> show the output ensemble mean and standard deviation estimates of the CAAE-ILUES and 
CAAE-DRDCN-ILUES methods for the 2-D and 3-D models, respectively.  
The reference output fields and the observation locations are also shown to facilitate the analysis of results. Notice that the statistics 
of the surrogate-based method are computed using the outputs predicted by the surrogate model rather than the forward model. 
It can be observed that the CAAE-DRDCN-ILUES method achieves similar ensemble mean estimates to those of the CAAE-ILUES method 
which successfully reproduce the main patterns of the reference output fields in the two cases.
Similar to the estimation of the conductivity field, a higher reproduction accuracy and a lower estimation uncertainty are observed near the observation wells than those in regions where no information is collected.  

The results indicate that the CAAE network is able to reconstruct well  non-Gaussian 
conductivity fields with different heterogeneity patterns and therefore the CAAE-ILUES method can obtain good inversion results in the 2-D and 3-D cases. However, for the two high-dimensional and highly-nonlinear 
inverse problems considered here, more than 60,000 forward model runs are needed in each case, leading to a high computational cost.  
On the contrary, the surrogate-based framework CAAE-DRDCN-ILUES can provide similar inversion results with much less computation costs. To assess the computational efficiency of the surrogate method, let N_C-I and N_C-D-I denote the number of forward model runs needed by CAAE-ILUES and CAAE-DRDCN-ILUES, respectively. Also, the C_f and C_train represent the computational cost of one forward model run and of training the surrogate model, respectively. The computation time of the original (C_C-I) and surrogate-based (C_C-D-I) ILUESs are written as 


    C_C-I=N_C-IC_f,


and


    C_C-D-I=N_C-D-IC_f+C_train,


respectively (the training time of CAAE parameterization model is the same for the two methods, thus it is not included in the equations). Thus, C_C-I=1.74×10^5 s and C_C-D-I=1.55×10^4 s in the 2-D case (N_C-I=62,000, N_C-D-I=4,000, C_f=2.8 s, and C_train=4.32×10^3 s), and C_C-I=8.25×10^5 s and C_C-D-I=9.02×10^4 s in the 3-D case (N_C-I=63,000, N_C-D-I=4,000, C_f=13.1 s, and C_train=3.78×10^4 s). The computational savings obtained by the surrogate-based method (C_C-I-C_C-D-I/C_C-I×100%) in the 2-D and 3-D cases are both about 90%. When applying the method to computationally more intensive forward models (i.e., when C_f is large), more savings of computation time 


    C_saving=C_C-I-C_C-D-I=(N_C-I-N_C-D-I)C_f-C_train,


can be expected as N_C-I≫ N_C-D-I.
 



§ CONCLUSIONS

In this study, we propose an integrated inversion framework for efficient characterization of solute transport in non-Gaussian conductivity fields with multimodal distributions. 
In the proposed framework, a CAAE network is developed for parameterization of the conducitivity field 
using a low-dimensional latent representation. In addition, a DRDCN network is developed for surrogate modeling 
of the solute transport model with high-dimensional and highly-complex input-output mappings. The two networks 
are combined with the ILUES inversion method <cit.> to formulate the CAAE-DRDCN-ILUES 
inversion framework. To improve the networks' performance for approximating the highly-complex mappings in the 
problems considered, in both network architectures, 
we adopt a multilevel residual learning structure referred to as residual-in-residual dense 
block <cit.>. 
This structure can effectively ease the training process allowing for a large network depth which has the potential to substantially increase the network's capability to approximate highly-complex mappings. 


The performance of the proposed method is evaluated using  2-D and 3-D solute transport modeling with
non-Gaussian conducitivity fields that have different patterns of conductivity heterogeneity. The results indicate that the CAAE is capable of representing a non-Gaussian conductivity 
field using a low-dimensional latent vector, though this is at the cost of the smoothness of local heterogeneity. The DRDCN network shows a superior performance over the 
DDCN network proposed in our previous studies <cit.> in obtaining accurate surrogate models.
The residual-in-residual dense block structure greatly improves the network's capacity in 
approximating highly-complex mappings. The application of the CAAE-DRDCN-ILUES method for estimation of 
the 2-D and 3-D non-Gaussian conductivity fields shows that it can obtain similar 
inversion results and predictive uncertainty estimations 
to those obtained by the original inverse method without surrogate modeling. The integrated 
method is highly efficient since the training of the surrogate model requires only a small number of forward model runs. 
The solute transport models considered in this work are relatively fast in order to quickly test the proposed method in a reasonable time. When applying 
to computationally more intensive forward models, significant computational gains can be expected. 

Note that incorporating the surrogate model in the inversion introduces an additional source of uncertainty due to approximation errors. Consideration of these approximation errors in the inversion process can be beneficial in improving the accuracy and reliability of the results <cit.>. To this end, even if not considered or demonstrated herein, one can build a  model for the approximation errors <cit.> or employ a Bayesian version of the DRDCN surrogate <cit.> for estimating the variance in the prediction, so as to incorporate the approximation errors in the inversion formula.
In the current work, the presented methods were demonstrated using synthetic problems. Their potential use in practical applications and other complex systems beyond groundwater solute transport deserves further exploration due to their data-driven and nonintrusive nature.

 

S.M. acknowledges partial financial support from the Center for 
Informatics and Computational Science (CICS) at the University of Notre Dame where this work was performed. The work of N.Z. was supported 
from the Defense Advanced Research Projects Agency under the Physics of Artificial Intelligence 
program (contract HR00111890034). The computing was supported via an AFOSR-DURIP award with additional   resources   provided by the University of Notre Dame Center for 
Research Computing. The work of X.S. and J.W. was supported by the 
National Natural Science Foundation of China (No. 41730856 and 41672229). The authors would like to thank the three anonymous reviewers and the Editors for their helpful comments.
The codes and data used in this work are available at <https://github.com/cics-nd/CAAE-DRDCN-inverse>.






§ DEEP DENSE CONVOLUTIONAL NETWORK

The deep dense convolutional network (DDCN) for surrogate modeling <cit.> 
is based on a dense connection structure called dense block <cit.>. The dense block 
introduces connections between its internal nonadjacent layers to enhance the information propagation 
through the network, so as to reduce the training sample size 
required to obtain desired approximation accuracy <cit.>. An illustration of the dense block 
structure is shown in Figure <ref>a. The difference between the dense blocks used 
in the DDCN network and in the DRDCN network proposed in this study is that, in DDCN, the input 
feature maps of the dense block's last layer are concatenated to its output feature maps to be 
fed into the next layer; while in DRDCN only the output feature maps of the dense block's last layer are 
passed to the next layer to allow an element-wise addition operation 
in the residual learning strategy (see section <ref>). 

We adopt the DDCN network architecture employed in our previous study <cit.> 
which was used for surrogate modeling of a 2-D solute transport model with Gaussian conductivity fields. The ReLU instead of Mish activation was used in DDCN.
The network is composed of 27 convolutional layers with three dense blocks. For the 3-D case considered in this study, we directly replace the 2-D convolutional layers in the network with the 3-D convolutional layers. 
More details about the network architecture can be found in mo2019inverse. When training the DDCN network, we use the same settings as in the DRDCN network. That is, the network is trained using the L_1 loss function defined in equation (<ref>) for 200 epochs  in the 2-D case and 300 epochs in the 3-D case  with the Adam optimizer <cit.>. The batch size is 32 and the initial learning rate is 5× 10^-3. A learning rate scheduler which drops ten times on plateau during training is used.  
Python codes of DDCN are available at <https://github.com/cics-nd/cnn-inversion>.

    


 

