Eric Laloy
email: elaloy@sckcen.be

Romain Hérault

John Lee

Diederik Jacques

Niklas Linde

Inversion using a new low-dimensional representation of complex binary geological media based on a deep neural network

published or not. The documents may come

Inversion using a new low-dimensional representation of complex binary geological media based on a deep neural network

Introduction

Inverse modeling plays a fundamental role in subsurface hydrology and many other Earth science disciplines. Basically, one iteratively proposes new model perturbations that are consistent with a given prior model until the resulting forward response agrees with a set of measured data up to a pre-specified level. The inversion outcome commonly consists of one or more 2D or 3D subsurface property field(s) describing, for instance, hydraulic conductivity or porosity. Due to conceptual and numerical errors in the forward problem formulation, (input and output data) measurement errors and insufficient information content in the measured data, the inverse problem may not admit a unique solution.

Hence, the inversion process should ideally provide an ensemble of geologic model realizations that accurately span the range of possible geologic models that honor both the available data used in the inversion and prior information.

The prior model is often conveniently specified by a two-point geostatistical model or variogram, by assuming the subsurface structure to be multi-Gaussian. Whenever this simplifying assumption does not hold (e.g., [START_REF] Gómez-Hernández | To be or not to be multi-Gaussian? A reflection on stochastic hydrogeology[END_REF][START_REF] Journel | The necessity of a multiple-point prior model[END_REF], the prior model can instead be informed by a so-called training image (TI). The TI is a large gridded 2D or 3D unconditional representation of the expected target spatial field that can be either continuous or categorical (e.g., geologic facies image). To generate model realizations that reproduce higher-order statistics found in the TI (together with possible direct point conditioning data), one can use a multiple-point statistics (MPS) simulation method (e.g., [START_REF] Guardiano | Multivariate geostatistics: beyond bivariate moments[END_REF][START_REF] Strebelle | Conditional simulation of complex geological structures using multiple point statistics[END_REF][START_REF] Hu | Multiple-point geostatistics for modeling subsurface heterogeneity: a comprehensive review[END_REF]Mariethoz et al., 2010b;[START_REF] Tahmasebi | Multiple-point geostatistical modeling based on the cross-correlation functions[END_REF].

Since MPS does generally not rely on a explicit mathematical description of the prior model in terms of model parameters, it is impossible to use it for classical parameterbased inversion. This prohibits the use of gradient-based methods such as randomized maximum likelihood (RML, [START_REF] Kitanidis | Quasi-linear geostatistical theory for inversing[END_REF][START_REF] Oliver | Conditioning permeability fields to pressure data[END_REF] or regular Markov chain Monte Carlo (MCMC) sampling [START_REF] Robert | Monte Carlo statistical methods[END_REF][START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF]. Note that this limitation also holds for the Markov mesh model (MMM, [START_REF] Stien | Facies modeling using a Markov mesh model specification[END_REF]) method despite the fact that it builds a statistical model describing transition probabilities of the facies distribution. Instead, one needs to resort to sampling from the prior model using MPS simulation. This process that is often referred to as sequential geostatistical resampling (SGR) inversion [START_REF] Hansen | Using geostatistics to describe complex a priori information for inverse problems[END_REF]Mariethoz et al., 2010a;[START_REF] Hansen | Inverse problems with non-trivial priors: efficient solution through sequential Gibbs sampling[END_REF] basically consists of iteratively generating new random perturbations from the prior model and to accept or reject them based on the resulting data misfit. The SGR inversion method is powerful in finding models that fit the data, but it suffers from two important drawbacks. As shown by [START_REF] Laloy | Merging parallel tempering with sequential geostatistical resampling for improved posterior exploration of high-dimensional subsurface categorical fields[END_REF], (1) it tends to only explore the immediate vicinity of a single global minimum even when the objective function/posterior landscape is highly multi-modal and (2) may produce degraded geologic model realizations when fitting large datasets with high signal-to-noise ratios (SNR).

Although some recent work has demonstrated substantial improvements to SGRbased inversion [START_REF] Zahner | Image synthesis with graph cuts: a fast model proposal mechanism in probabilistic inversion[END_REF], in this study we take another direction and propose a new parametric dimensionality reduction approach for complex binary prior models. The driving idea is that if one can build a lower-dimensional model parameterization from which one can sample (after appropriate non-linear transformations) model realizations that are consistent with the TI, then global parameter-based inverse methods for moderately large continuous parameter spaces can be used. For instance, this opens up the possibility to explore the posterior model distribution using state-of-the-art adaptive MCMC sampling with DREAM (ZS) [START_REF] Vrugt | Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling[END_REF]Laloy and Vrugt , 2012).

Previous work on using compressed parametric bases for hydrogeological inversion has relied on principal component analysis (PCA, e.g., [START_REF] Reynolds | Reparameterization techniques for generating reservoir description conditioned to variograms and well-test pressure data[END_REF][START_REF] Sarma | Efficient real-time reservoir management using adjoint-based optimal control and model updating[END_REF], kernel-PCA (e.g., [START_REF] Sarma | Kernel principal component analysis for efficient, differentiable parameterization of multipoint geostatistics[END_REF], level-set (e.g., [START_REF] Dorn | History matching of petroleum reservoirs using a level set technique[END_REF], discrete wavelet transforms (e.g., [START_REF] Awotunde | Reservoir description with integrated multiwell data using two-dimensional wavelets[END_REF], discrete cosine transform (DCT, e.g., [START_REF] Jafarpour | Compressed history matching: exploiting transform-domain sparsity for regularization of nonlinear dynamic data integration problems[END_REF][START_REF] Linde | Distributed soil moisture from crosshole ground-penetrating radar travel times using stochastic inversion[END_REF], singular value decomposition (SVD, e.g., [START_REF] Tavakoli | Monte Carlo simulation of permeability fields and reservoir performance predictions with SVD parameterization in RML compared with EnKF[END_REF], and K-SVD (e.g., [START_REF] Khaninezhad | Sparse geologic dictionaries for subsurface flow model calibration: part I. Inverse formulation[END_REF][START_REF] Khaninezhad | Sparse randomized maximum likelihood (SpRML) for subsurface flow model calibration and uncertainty quantification[END_REF] to name most of the strategies used. However, none of these bases define a manifold that is restricted to complex (non-Gaussian) geological models that are in strong agreement with a specific TI. In other words, randomly sampling the parameter space produces model realizations that are generally inconsistent with the TI. To the best of our knowledge, there is so far only the optimization-PCA approach (O-PCA) by Vo andDurlofsky (2014, 2015) that aims to build a lower-dimensional model representation that is fully consistent with the TI. The O-PCA essentially consists of a post-processing of a PCA model. To work it is necessary that the underlying PCA representation encodes sufficient information regarding important patterns in the TI model. As illustrated later on in this paper, this may not be the case when the TI contains highly connected features with a relatively large degree of variability.

In this work, we use for the first time a deep neural network (DNN, see, e.g., [START_REF] Goodfellow | Deep learning[END_REF], for an overview) to build a parametric low-dimensional representation of complex geologic models. More specifically, we train the generator of a deep variational autoencoder (VAE, [START_REF] Kingma | Auto-encoding variational Bayes[END_REF] such that it generates TI-consistent geologic model realizations when fed with a low-dimensional standard Gaussian noise vector.

A key characteristic of our VAE is that it is largely made up of convolutional layers (see, e.g., [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF][START_REF] Krizhevsky | ImageNet classification with deep convolutional neural networks[END_REF]. This type of neural layer is based on an adaptive convolutional filter and is well suited to image processing [START_REF] Krizhevsky | ImageNet classification with deep convolutional neural networks[END_REF][START_REF] Goodfellow | Deep learning[END_REF]. Our model generator is then used within a Bayesian inversion framework to sample the posterior distribution of binary 2D and 3D subsurface property fields.

This paper is organized as follows. Section 2 presents the different elements of our dimensionality reduction (DR) approach and the selected inversion framework. This is followed in section 3 with the analysis of its performance in generating geostatistical realizations using 2D and 3D TIs with and without conditioning on direct point data (also referred to as hard data), together with some comparisons against other existing DR methods. Synthetic 2D and 3D experiments involving both steady-state and transient groundwater flow are then used to demonstrate our proposed DR-based inversion approach. In section 5 we discuss the advantages and limitations of our method and outline possible future developments. Finally, section 6 concludes with a summary of the most important findings.

Methods

Deep neural network architecture

Generalities

We consider a deep neural network that belongs to the class of autoencoders (AEs, see, e.g., [START_REF] Goodfellow | Deep learning[END_REF]. Neural networks basically define the (possibly complex) relationships existing between input, x, and output, y, data vectors by using combinations of computational units that are called neurons. A neuron is an operator of the form:

h (x) = f ( x, w + b) , (1)
where h (•) is the scalar output of the neuron, f (•) is a nonlinear function that is called the "activation function", •, • signifies the scalar product, w = [w 1 , • • • , w N ] is a set of weights of same dimension, N , as x and b represents the bias associated with the neuron. For a given task, the values for w and b associated with each neuron must be optimized or "learned" such that the resulting neural network performs as well as possible. When f (•) is differentiable, w and b can be learned by gradient descent. Common forms of f (•) include the rectified linear unit (ReLU), sigmoid function and hyperbolic tangent function.

When there is no directed loops or cycles across neurons or combinations thereof, the network is said to be feedforward (FFN). In the FFN architecture, the neurons are organized in layers. A standard layer is given by

h (x) = f (Wx + b) , (2)
where W and b are now a matrix of weights and a vector of biases, respectively. The name multilayer perceptron (MLP) designates a FFN with more than one layer. A most typical network is the 2-layer MLP, which consists of two layers with the outputs of the first-layer neurons becoming inputs to the second-layer neurons

y = g [h (x)] ≡ f 2 [W 2 • f 1 (W 1 x + b 1 ) + b 2 ] , (3)
where g (•) and h (•) are referred to as output layer and hidden layer, respectively.

In theory, the two-layer MLP described in equation ( 3) is a universal approximator as it can approximate any underlying process between y and x [START_REF] Cybenko | Approximations by superpositions of sigmoidal functions[END_REF][START_REF] Hornik | Approximation capabilities of multilayer feedforward networks[END_REF]. However, this only works if the dimension of h (•) is (potentially many orders of magnitudes) larger than that of the input x, thereby making learning practically infeasible and the two-layer MLP approximator useless for large N (typically N ≥ 10-25). For high-dimensional input data such as images, researchers have found that it is much more efficient to use many hidden layers rather than increasing the size of a single hidden layer (e.g., [START_REF] Goodfellow | Deep learning[END_REF]. When a FFN/MLP has more than one hidden layer it is considered to be deep. Nevertheless, current deep networks are not necessarily purely FFN but may mix different aspects of FFN, such as convolutional neural networks (CNN, see section 2.1.3 below) and recurrent neural networks (RNN, see, e.g., [START_REF] Goodfellow | Deep learning[END_REF].

Variational autoencoders

An AE is a deep neural network that defines a reversible, nonlinear low-dimensional parameterization of (higher dimensional) input data. Consequently, it has an hourglasslike shape in terms of neural layer sizes (see simplified representation in Figure 1). The most central layer is referred to as the "code" and defines the low-dimensional space, while the output layer has the same dimensionality as the input layer. The part of the network that connects the input to the code is called the "encoder", while the one that connects the code to the output is referred to as the "decoder".

Autoencoders are generative, which means that, after appropriate training (see section 2.2), they can be used to generate new pattern realizations that are consistent with those found in a given set of (training) features. However, only the class of AEs formed by the so-called variational autoencoders (VAE) can use white noise as input to generate new patterns. In a geostatistical context, this implies that the VAE can be trained such that it randomly generates new geologic model realizations that honor the higher-order statistics found in a set of training images. How to achieve this is further described later on in this section and in section 2.2.

A VAE network can be summarized as follows. The code consists of two lowerdimensional vectors: a vector of means, µ, and a vector of standard deviations, σ (Figure 2). The first element of the decoder is a randomly sampled vector of "latent" standard normal variables, z l , of the same size, d, as µ and σ. The z l vector is subsequently rescaled into z = z l × σ + µ (where × means element-wise multiplication), and the soproduced z vector continues its journey through the decoder to eventually produce the output, x (Figure 2). Here the symbol x is used instead of y to make it clear that for this application the output is a reconstruction of the input, x. As detailed in section 2.2, the w vectors and b values that are associated with every neuron of every layer of the network (encoder part, µ, σ and decoder part) are jointly optimized by gradient descent such that (1) the differences between x and x are minimized and (2) the d-dimensional z vector conforms as closely as possible to a N (0 d , I d ) with 0 d a d-dimensional zero-vector and I d the d × d identity matrix. After training, new model realizations can be generated by sampling z ∝ N (0 d , I d ) and running the z samples through the decoder.

At this point, we would like to re-emphasize that it is the combination of (1) a lowdimensional parameterization in terms of a multivariate standard normal distribution and (2) the ability to map a given set of low-dimensional parameter values into complex MPS geostatistical model realizations that makes our presented method unique and suitable for MPS-based inversion. Figure 1: Schematic example of an autoencoder with three hidden layers, h 1 , h 2 and h 3 . The x = [x 1 , • • • , x 6 ] vector is a 6-dimensional input feature, the h (•,•) denote the hidden layers with (red-contoured) h (2,•) the central layer, or "code", where the largest dimensionality reduction is performed and the x = [x 1 , • • • , x6 ] vector signifies the reconstructed input.

Figure 2: Illustration of a variational deep autoencoder (VAE). The boxes labeled with the µ and σ letters represent the low-dimensional mean and standard deviation vectors, respectively, that are encoded from the input data. During training, z l is drawn from the fixed p (z l ) distribution and then rescaled into z = z l × σ + µ (where × means elementwise multiplication). Training aims at jointly minimizing (1) the data reconstruction loss and the deviations between (2) µ and the zero-vector of the same size and (3) σ and the unit vector of the same size. The ensemble of tasks (2) and ( 3) is called the "reparameterization trick" [START_REF] Kingma | Auto-encoding variational Bayes[END_REF]. For subsequent stochastic data generation, the process starts with sampling z from the fixed p (z) distribution.

Convolutional layers

The last salient characteristic of our deep VAE is the use of convolutional layers at various levels of both the encoder and the decoder. The convolutional layer is the main building block of the CNN-type of architecture. It is particularly well suited for image processing applications because it explicitly accounts for the spatial structure of the input data, whether the input image, X, or the incoming hidden layer in the network. When the input is a 2D image (with possibly 3 channels for a RGB image), a convolutional layer,

h k , is built from a series of k = 1, • • • , N k small N i × N j filters, w k , that convolve an input pixel, X m,n to h k m,n as h k m,n (X m,n ) = f   N i ,N j i=1,j=1 w k i,j X m+i,n+j + b k   , (4)
where for computational efficiency f (•) is typically a ReLU: f (x) = max (0, x). The so-produced ensemble of N k h k "feature maps" forms a volume called the convolutional representation. The larger the N k , the potentially more comprehensive is the representation of the input data. Other important convolution parameters that we do not discuss here are the "stride", that is, the degree of overlapping between successive moves in the forward pass of a given filter, and the "zero-padding" which involves padding the borders of the input image or volume for size preservation. For further information, we refer the reader to [START_REF] Goodfellow | Deep learning[END_REF] and online tutorials1 .

Note that each neuron of a given convolutional layer (volume) is connected to only a local region of the input image or volume, called "receptive field" (see Figure 3). This enforces spatially-local exploration at the level of one convolutional layer. The encoded spatial information includes increasingly larger patterns when the number of stacked convolutional layers increases. Also, recall that the elements (pixels/voxels) of a given feature map share the same w k and b k . The number of CNN parameters to be optimized is thus greatly reduced compared to that of a standard fully-connected (FC) (or dense) FFN of similar depth. Besides the convolutional layer itself, the other important component of a CNN is the "pooling" layer. Pooling is a form of nonlinear down-sampling used to reduce the dimensions of a convolutional layer and thereby limit overfitting. It proceeds by converting small sub-regions of equal size into single values, such as the maximum or the mean. Using the maximum ("max-pooling") is the most common practice. When used, pooling always follows a series of one or more convolution operations.

Figure 3: Stylized representation of convolutional layer neurons connected to their receptive field (picture inspired by http://cs231n.github.io/convolutional-networks/).

Practical implementation

We now have described all of the main ingredients of our devised deep VAECNN network, which is depicted in Figure 4. For the sake of brevity, the various individual fullyconnected (FC), variational (Q), convolutional (C) and pooling (P) layers (see Figure 4) are not detailed in this paper but the corresponding computer code is available from the first author. Our VAECNN was implemented within the open-source LASAGNE Python software (Lasagne, 2017) which works on top of the open-source THEANO Python library [START_REF] Development | THEANO: A Python framework for fast computation of mathematical expressions[END_REF]. In addition, note that for each result presented in this study, the z vector is always 50-dimensional no matter the dimensions of the considered geologic model domain (from 100 × 100 to 30 × 32 × 27).

Unless stated otherwise, model generation is performed by recycling the initially produced model (bottom blue architecture in Figure 4) through the full network (top light red architecture in Figure 4) for 10 times. This sequential relooping generally improves the quality of the generated model realizations. The very last step of our generation process is an hard thresholding of the produced model. We used the 0.5 cutoff value although the classical Otsu's method could be used as well [START_REF] Otsu | A threshold selection method from gray-Level histograms[END_REF]. Fortunately, the impact of relooping on generation time remains quite acceptable. For the models considered herein, generating a realization without relooping takes about 0.01 s -0.03 s on a modern workstation. Including the 10 successive cycles then incurs a total generation time in the range 0.2 s -0.4 s.

Training the deep variational autoencoder

Considering the case where a 100 × 100 binary image is generated from a 50-dimensional z, the network described in Figure 4 includes 34,039,237 optimized parameters (weights and bias associated with neurons of dense layers and filters of convolutional layers). Such a high number of optimized parameters is not unusual for applications of deep neural networks. As the activation functions are differentiable, using an objective function that is also differentiable makes it possible to optimize/train the network by stochastic gradient descent (that is, gradient descent using a series of mini-batches rather than all the data at once) together with back propagation. This means that the objective function derivative is propagated backwards throughout the network using the chain rule, in order to update the parameters. Various stochastic gradient descent algorithms are available. In this work, we used the adaptive moment estimation (ADAM) algorithm which has been proven efficient for different types of deep networks [START_REF] Kingma | Adam: a Method for Stochastic Optimization[END_REF].

The objective or loss function, L, must be set such that at training time: (1) the differences between reconstructed images, X, and their original counterparts, X, are minimized and (2) the d-dimensional random z vector conforms as closely as possible to a N (0 d , I d ). This translates into the following formulation [START_REF] Kingma | Auto-encoding variational Bayes[END_REF][START_REF] Gregor | DRAW: A recurrent neural network for image generation[END_REF])

L = L x + L z , (5)
with

L x = -log [p (X|z)] , (6)
L z = D KL [q (z|X) ||p (z)] , (7)
where p (z) is the probability density function (pdf) of the code z (here p (z) ≡ N (0 d , I d )) and the conditional distribution p (X|z) is referred to as a probabilistic decoder, since given z it generates a probability distribution over the possible corresponding values of X. Similarly, q (z|X) denotes a probabilistic encoder, as given a feature X it creates a distribution over the possible values of z from which X could have been produced [START_REF] Kingma | Auto-encoding variational Bayes[END_REF]. To ensure that q (z|X) is close to the targeted p (z), the Kullback-Leibler divergence (D KL ) from q (z|X) to p (z), D KL [q (z|X) ||p (z)] is minimized (equation ( 7)).

If X is a binary image, equation ( 6) can be rewritten as a binary cross-entropy minimization

L x = N i=1 -X i log Xi -(1 -X i ) log 1 -Xi . (8)
In addition, for p (z) ≡ N (0 d , I d ) equation 5 becomes

L z = 1 2 d i=1 µ 2 i + σ 2 i -log σ 2 i - d 2 , ( 9
)
where the µ i and σ i correspond to the i th elements of the µ and σ vectors displayed in Figure 2, respectively. Working with µ and σ vectors rather than with the z vector directly is called the "reparameterization trick" (for details, see [START_REF] Kingma | Auto-encoding variational Bayes[END_REF].

Finally, in practice one may need to weight the components of L, leading to

L = L x + αL z , ( 10
)
where α is a weight factor that we set after limited trial and error either to 20 or to 40 depending on the application.

The minimization of equation ( 10) was performed on a GPU Tesla K40 and training the VAECNN for 100 epochs (full cycles of the stochastic gradient descent) took between 5 and 13 hours, depending on the size of the model domain (from 100×100 to 32×30×27), the number of training images (from 19,000 to 80,000), and the exact size of the network. Computing times will be further detailed below for each test case individually.

Assessing geostatistical simulation quality

Although inversion is our primary objective, assessing the quality of the geostatistial realizations produced by our DR approach is important as it will determine the ultimate quality of our inversion results. Following [START_REF] Tan | Comparing training-image based algorithms using an analysis of distance[END_REF], we considered both how well the patterns in the training images are reproduced by the DR-based realizations and the between-realization variability, also termed "space of uncertainty".

Besides visual inspection, we estimated closeness with the training images by means of the two-point cluster function (CF, Torquato et al., 1988) which is also called connectivity function [START_REF] Pardo-Igúzquiza | CONNEC3D: A computer program for connectivity analysis of 3D random set models[END_REF]. The CF is the probability that there exists a continuous path of the same facies between two points of the same facies separated by a given lag distance (see, [START_REF] Torquato | Two-point cluster function for continuum percolation[END_REF][START_REF] Pardo-Igúzquiza | CONNEC3D: A computer program for connectivity analysis of 3D random set models[END_REF], for mathematical details). Using the implementation by Lemmens et al. (2017), the CF was calculated for each facies along the x, y, and main diagonal, d xy , directions for the 2D case studies. For the 3D case studies, we also considered the z and main diagonal d xz and d yz directions. This resulted into a total of six CF curves for each combination of indicator and facies.

The so-called space of uncertainty is a measure of the between-realizations variability associated with given geostatistical simulation algorithm. We characterized the space of uncertainty for the 2D case using an average distance between multiple-point histograms (MPH) over 100 realizations (see Tan et al., 2014, for details). This was done both for the training set and DR-based realizations. The used distance is the so-called Jensen-Shannon (JS) divergence [START_REF] Cover | Elements of information theory[END_REF][START_REF] Tan | Comparing training-image based algorithms using an analysis of distance[END_REF].

d JS m k , m k , = 1 2 N b i=1 M P H k i log M P H k i M P H k i + 1 2 N b i=1 M P H k i log M P H k i M P H k i , (11)
where m k and m k are two different model realizations and

M P H k 1 , • • • , M P H k N b
are the components of the MPH k vector calculated for realization m k . Since the MPH was computed for a 4 × 4 template, N b was 2 16 = 65, 536.

The space of uncertainty, defined here as the average distance, d JS , for the set of either K training images or K DR-based realizations, is given by

d JS = 1 K (K -1) K k=1 K k =1 d JS m k , m k , . (12)

Bayesian inversion

A common representation of the forward problem is

d = F (θ) + e, (13)
where

d = (d 1 , . . . , d N ) ∈ R N , N ≥ 1 is the measurement data, F (θ)
is a deterministic forward model with parameters θ and the noise term e lumps all sources of errors.

In the Bayesian paradigm, parameters in θ are viewed as random variables with a posterior pdf, p (θ|d), given by

p (θ|d) = p (θ) p (d|θ) p (d) ∝ p (θ) L (θ|d) , (14)
where L (θ|d) ≡ p (d|θ) signifies the likelihood function of θ. The normalization factor p (d) = p (θ) p (d|θ) dθ is not required for parameter inference when the parameter dimensionality is fixed. In the remainder of this paper, we will thus focus on the unnormalized density p (θ|d) ∝ p (θ) L (θ|d).

To avoid numerical over-or underflow, it is convenient to work with the logarithm of L (θ|d) (log-likelihood): (θ|d). If we assume e to be normally distributed, uncorrelated and with known constant variance, σ 2 e , (θ|d) can be written as

(θ|d) = - N 2 log (2π) -N log (σ e ) - 1 2 σ -2 e N i=1 [d i -F i (θ)] 2 , ( 15
)
where the F i (θ) are the simulated responses that are compared the

i = 1, • • • , N measurement data, d i .
An exact analytical solution of p (θ|d) is not available for the type of non-linear inverse problems considered herein. We therefore resort to MCMC simulation (see, e.g., [START_REF] Robert | Monte Carlo statistical methods[END_REF]. More specifically, the DREAM (ZS) algorithm is used to approximate the posterior distribution. A detailed description of this sampling scheme including a proof of ergodicity and detailed balance can be found in [START_REF] Vrugt | Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling[END_REF] and Laloy and Vrugt (2012). Multiple contributions in hydrology and geophysics (amongst others) have demonstrated the ability of DREAM (ZS) to sample target distributions with 50-250 dimensions (Laloy et al., 2012[START_REF] Laloy | Efficient posterior exploration of a high-dimensional groundwater model from two-stage Markov Chain Monte Carlo simulation and polynomial chaos expansion[END_REF][START_REF] Linde | Distributed soil moisture from crosshole ground-penetrating radar travel times using stochastic inversion[END_REF][START_REF] Laloy | Probabilistic inference of multi-Gaussian fields from indirect hydrological data using circulant embedding and dimensionality reduction[END_REF][START_REF] Lochbühler | Summary statistics from training images as prior information in probabilistic inversion[END_REF].

Geostatistical simulation results

Before considering our main goal, that is, solving inverse problems (section 4), we investigate the performance of our DR approach for unconditional and conditional geostatistical simulation. This step is important as we aim at producing inversion models that not only honor the data used, but also the spatial statistics of our selected TI. One thus need to make sure that model realizations produced by randomly sampling the constructed low-dimensional base are consistent with the TI.

2D channelized aquifer

Our first example considers unconditional model generation of a 100 × 100 2D binary and channelized aquifer. The underlying TI is closely related to the classical 250 × 250 binary training image introduced by [START_REF] Strebelle | Conditional simulation of complex geological structures using multiple point statistics[END_REF] that is routinely used to test new MPS algorithms. Our used TI has some slight additional variability in that the channel rotation angles are allowed to deviate (at most) from -15°to +15°compared to the rotation angles found in the original TI. The training set used to train the VAECNN contains 80,000 100 × 100 geostatistical realizations generated with the direct sampling (DS) approach as implemented in the DeeSse MPS algorithm (Mariethoz et al., 2010b). Obtaining these 80,000 realizations took 14 hours using eight CPUs in parallel. The choice of 80,000 realizations was dictated by available computational resources. Limited testing further indicated that using 40,000 model realizations for training would very likely not change the results presented below. Training the VAECNN for this example incurred a computational time of approximately 13 hours. As z is 50-dimensional, the achieved compression ratio is 200.

Figure 5 presents four (randomly chosen) MPS realizations from the training set, together with eight (randomly chosen) model realizations obtained by our DR approach. The DR-based realizations resemble the training model realizations well, although some artifacts remain in terms of a moderately larger occurrence of broken channels and a slight oversmoothing of the channels. Figure 6 Varying this number in the range 50 -10000 did not improve generation performance. For DCT, 250 coefficients were used as using a smaller number of coefficients leads to overly degraded models in direct compression mode (compression of existing models). The DCT-based model realizations were then generated by randomly sampling uniform distributions defined by the empirical upper and lower bounds of the 250 maximum (in absolute value) coefficients associated with the training set. The DCT realizations were thresholded such that the facies proportions are close to that of the TI (facies 0: 0.7, facies 1: 0.3). Varying the number of DCT coefficients has no significant impact on generation quality. Clearly, none of these three approaches (PCA, OPCA, DCT) sample from a proper low-dimensional parameter space. In contrast, our proposed approach does a much better job. Note that PCA, OPCA and DCT are far from providing stateof-the-art MPS simulations. Instead, up to now they formed the state-of-the-art in terms of parametric low-dimensional model representations that have been used with advanced inversion methods. The key of our method a low-dimensional parameterization that is able to provide high-quality MPS realizations (compare Figures 5 and7). Our second example considers a case where direct conditioning data are available. This example is similar to example 1 (Figure 5) except that facies at 9 locations are known (red circles in Figure 8). To handle direct conditioning, our approach requires that the MPS-generated realizations forming the training set honor the conditioning points. Building a training set of 40,000 conditioned models with DS took 7-8 hours using again eight CPUs in parallel. The training of our VAECNN then took 6.5 hours.

A set of (randomly chosen) training and DR-generated model realizations are shown in Figure 8. A similar generation performance is visually observed as for example 1 (Figure 5) and the DR-based realizations show again similar CF curves as the training images (Figure 9). However, the achieved space of uncertainty is 12% smaller than that of the conditioned training set:

d JS (DR) d JS (T R) = 0.88.
The conditioning accuracy is assessed by analyzing 1000 random model realizations. We find that 68% of the model realizations honor all of the 9 prescribed facies and that 97% of the model realizations contain at most one mismatching conditioning datum. Also, the more common facies 0 is more frequently honored (98%) than the less frequent facies 1 (91%). This is less than for our DS-generated training set that always honor the nine conditioning points. Nevertheless, the training realizations (Figures 8a-d) often contain 1 or 2 isolated conditioning point(s) within an homogeneous zone of the other facies. This type of artifact is less frequent in the DR-based model realizations (Figures 8e-l). The main reason for the conditioning errors is the sequential relooping used for producing a realization (see section 2.1.4). Without relooping, the model realizations present more frequently broken channels and/or small groups of isolated pixels (not shown), but the percentage of realizations that fully honor the 9 conditioning data is 92% (instead of 68% with relooping) and 100% (instead of 97% with relooping) of the realizations has at most a single erroneous datum.

3D model

We now turn our attention to the generation of 3D model realizations. The selected TI is the Maules Creek valley alluvial aquifer available from http://www.trainingimages. org/training-images-library.html. For the unconditional case, a training set of 19,500 30 × 32 × 27 model realizations was built using the recent graph cuts (GC) patch-based MPS algorithm by [START_REF] Li | Patch-based iterative conditional geostatistical simulation using graph cuts[END_REF]. Running 8 different GC instances simultaneously, this took approximately 6 hours. The GC method was selected because it is much faster than DS. In addition, training our VAECNN took about 7 hours. The convolutional layers in our VAECNN are designed for 2D images only, but with possibly different color channels. To generate 3D images, one can use as many channels as there are horizontal layers in the model ( 27). This does not explicitly account for patterns in the vertical direction since the convolutional operations are performed in the horizontal plane only. Nevertheless, equation ( 8) penalizes reconstructed models that deviate from the training set in all three spatial directions.

Figure 10 compares (randomly chosen) training and DR-based realizations. As the model domain contains 30×32×27 = 25, 920 voxels, the compression ratio is 25920/50 ≈ 518. The proposed DR approach creates realizations that are similar to those in the training set, although isolated voxels of each facies are over-represented. The realizations also reproduce well the CF metrics of the training set (Figure 11).

We also considered the same 3D example with 56 direct conditioning points. These points correspond to the locations of the screens along the vertical multi-level piezometers used in our synthetic 3D inverse problem described in section 4.2 (Blue line segments in Figure 12). A GC-based training set containing 19,000 conditioned model realizations was used for training our VAECNN. Running again 8 different GC instances simultaneously, building the training set took approximately 14 hours. Similarly as for the unconditional case, training then lasted for 7 hours. Using default algorithmic settings, these original GC simulations jointly honor all of the 56 conditioning data only 6% of the times but contain at most 10 mismatching point 94% of the times. The DR-based realizations are again visually close to the training set (not shown), with a much smaller tendency to overproduce isolated voxels of each facies than for the unconditional case (not shown) and CF statistics that closely match those of the training set (not shown). Furthermore, our DR approach is found to condition at a slightly higher level than the selected GC algorithm. Among 1000 random DR-based realizations, 7% of the realizations honor all of the 56 prescribed facies while 98% of the realizations contain a most 10 mismatching points. We attribute these slightly superior statistics compared to the training set to a combination of randomness and complex bias in the dimensionality reduction. Overall, our VAECNN is found to condition equally well as the GC algorithm used to produce the training set. Our first inversion case study considers steady-state flow within a channelized aquifer. The 100 × 100 aquifer domain lies in the xy plane with a grid cell size of 1 m and a thickness of 1 m. Channel material and matrix material (see Figures 14a and15a) are assigned hydraulic conductivity values of 1 × 10 -2 m/s and 1 × 10 -4 m/s, respectively. Steady state groundwater flow is simulated using MODFLOW 2005 [START_REF] Harbaugh | Geological Survey modular groundwater model -the groundwater flow process[END_REF] assuming no flow boundaries at the upper and lower sides and fixed head boundaries on the left and right sides of the domain to ensure a lateral head gradient of 0.01 (-), with water flowing in the x-direction. A pumping well extracting 0.001 m 3 /s is located at the center of the domain. Simulated heads are collected at 49 locations that are regularly spread over the domain (Figures 14a and15a). These data were then corrupted with a Gaussian white noise using a standard deviation of 0.02 m. For the selected white noise realization, the measurement data have a root-mean-square-error (RMSE) of 0.0200 m. The corresponding SNR, defined as the ratio of the average RMSE when drawing prior realizations with our DR algorithm to the noise level is 25. No direct conditioning (facies) data are used.

We first verify whether a model realization belonging to the 50-dimensional manifold can be retrieved by the inversion process. To do so, the reference model (Figure 14a) was randomly generated using our DR approach. The DREAM (ZS) was ran in parallel, using 8 interacting Markov chains distributed over 8 CPUs. Uniform priors in [-5, 5] were selected for the 50 DR variables, θ. This might seem to contradict the theory behind our VAECNN which is based on calibration to a standard normal distribution for p (z). However, using a standard normal prior for θ was found to restrict the sampled model space too much for the inversion process to produce model realizations that in average fit the data to prescribed noise level. Our explanation for this phenomena is that with 50 DR variables the MCMC simulation needs sufficient freedom in the explored model space to eventually produce consistent models. Also, and equally important, is the observation that after training using equations (8-9), the quality of the generated realizations by our VAECNN is similar when using a zero-mean uniform or a standard normal distribution for p (z). From a pragmatic point of view, the zero-mean uniform prior can thus be used.

The chains start to jointly sample the posterior distribution, p (θ|d), after a (serial) total of 120,000 iterations, that is, 15,000 parallel iterations per chain (Figure 13a). The sampled realizations closely resemble the true model and the posterior variability is overall small (Figure 14). To assess sampling accuracy, we calculate over the last 10,000 posterior realizations in each chain the average fraction of pixels with facies identical to that of the true model, f PO , which is 0.86. The same quantity for the prior, f PR , can be directly calculated as 0.7 × 0.75 + 0.3 × 0.25 = 0.60 with 0.7 and 0.3 the prior fractions of each facies in the TI and 0.75 and 0.25 the fractions of each facies in the true model. The f PO f PR ratio will be used to compare the examples considered in this section. Here

f PO f PR = 0.86 0.60 = 1.43.
After a total of 200,000 MCMC iterations, the [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF] convergence diagnostic, R, is satisfied (i.e., R ≤ 1.2) for 37 out of the 50 sampled parameters (see, e.g., [START_REF] Laloy | Probabilistic inference of multi-Gaussian fields from indirect hydrological data using circulant embedding and dimensionality reduction[END_REF], for details about the use of R with DREAM (ZS) ).

The MCMC sampling should therefore be continued for a longer time for (official) full exploration of the posterior distribution.

Our second and more comprehensive test uses a true model that was generated by the DS algorithm. If the DR-representation is inappropriate, this true model might not be part of the 50-dimensional DR space. Identical inversion settings were used as in the previous example. Here the 8 chains start to jointly sample p (θ|d) after some 12,500 iterations per chain (Figure 13b). The posterior model realizations look visually close to the reference model, though posterior variability appears to be a bit larger than for the previous case (Figure 15). The associated f PO and f PR values are 0.78 and 0.7 × 0.71 + 0.3 × 0.29 = 0.59, respectively, which gives a f PO f PR ratio of 1.32. The results displayed in Figure 15 inspire confidence that the proposed DR-based inversion approach can retrieve consistent posterior models. With respect to posterior exploration, formal

convergence is yet to be declared after a total of 200,000 serial iterations, with the R criterion being satisfied for 29 out of the 50 dimensions of θ. Full posterior exploration therefore requires longer chains. To benchmark our approach against state-of-the-art MPS inversion, we performed eight independent SGR runs in parallel for the same case study as above and a similar CPU budget of 25,000 MCMC iterations per run. Moreover, the SGR algorithmic settings are optimal (for the considered type of application) as defined in [START_REF] Laloy | Merging parallel tempering with sequential geostatistical resampling for improved posterior exploration of high-dimensional subsurface categorical fields[END_REF]. It is observed that 5 SGR trials out of 8 fit the data to the appropriate noise level after less than 10,000 iterations (not shown). Three trials even reach the targeted 0.02 m misfit after less than 3000 iterations (not shown). Yet three other SGR trials do not sample the posterior distribution after the 25,000 iterations (trials #3, #4 and #8 in Figure 16). Furthermore, all SGR trials end up in a relatively narrow local optimum (not shown) while none of them find a model that is visually similar to the true model (Figure 16). Many of the realizations sampled by SGR are also degraded compared to the prior model (see top row of Figure 5), with isolated patches and broken channels (Figure 16). This translates into a quite lower f PO f PR ratio than for our DR-based approach: 1.08 (the average f PO over the 8 SGR trials is 0.64). We thus conclude that for the considered example, our DR approach is a superior alternative to SGR.

Case study 2: 3D transient hydraulic tomography

Our last case study focuses on 3D transient hydraulic tomography (e.g., [START_REF] Cardiff | Hydraulic conductivity imaging from 3-D transient hydraulic tomography at several pumping/observation densities[END_REF]. Transient head variations caused by discrete multilevel pumping tests were simulated using MODFLOW2005 within a 3D confined aquifer of size 30 × 32 × 27 with a voxel size of 1 m × 1 m × 1 m. This aquifer consists of a part of the Maules Creek TI discussed in section 3.2. The true model is depicted in Figure 17a. The multilevel discrete pumping setup consists of a 27-m deep, central multi-level well in which water is sequentially extracted every 4 m along a 1-m long screen (at depths of 3, 7, 11, 15, 19, 23 and 27 m respectively) during 30 minutes at a rate of 20 liters/min. The locations of this multilevel pumping well and those of the 8 surrounding multilevel observation wells are displayed in Figure 12. For each pumping sequence, drawdown are recorded every 4 m along a 1-m long screen in the 8 multi-level piezometers (Figure 12). For each drawdown curve, data acquired at the same four measurement times were retained leading to a total of 8 × 7 × 7 × 4 = 1568 measurement data. These measurement times were considered to be the four most informative ones after visual inspection of several drawdown curves (not shown). These data were corrupted with a Gaussian white noise using again a standard deviation of 0.02 m, which induced for the selected white noise realization a RMSE of 0.0200 m. The associated SNR is 6. In addition, the facies of the 56 locations (voxels) where pressure head is measured (blue line segments in Figure 12) are known and used as direct conditioning data in the inversion.

For this case study, the DREAM (ZS) sampler evolves 16 Markov chains in parallel using 16 CPUs. Uniform priors in [-5, 5] were again selected for the 50 dimensions of θ. Within the allowed computational budget of 60,000 iterations per chain, the 16 chains converge towards a data misfit in the range of 0.0202 m -0.0204 m (Figure 13c). This interval is close to the target level of 0.0200 m. This indicates that the posterior mode has not been sampled yet. The variability among the sampled models is rather small (Figure 17). The latter is likely due to the combination of two factors: (1) the peakedness of the likelihood function (equation ( 15)) caused by the large number of measurement data, and (2) the difficulties encountered by the sampling algorithm for exploring this complex target distribution. Notwithstanding, the sampled models (Figure 17b-i) are visually close to the true model (Figure 17a).

Discussion

Our proposed DR approach that is based on a deep neural network works well for the considered inversion of 2D steady-state flow data. It also provides useful solutions to our inverse problem involving 3D transient flow. However, further improvements are needed. As for most compression approaches, our dimensionality reduction scheme builds a lowdimensional representation where the value in each dimension influences every pixel/voxel of the generated realization in the original space. This "non-localization" poses problems for inverse problems involving categorical fields. Indeed, [START_REF] Laloy | Merging parallel tempering with sequential geostatistical resampling for improved posterior exploration of high-dimensional subsurface categorical fields[END_REF] have shown for SGR that iterative perturbation of a small well defined fraction of the current model works much better than slightly perturbing the whole model at every step. The same logic applies to compression-based inversion. Indeed, Figure 18 demonstrates for our first inversion case study that (1) changing the value of a single low-dimensional parameter influences the whole model realization domain and (2) even a small variation in the model realization can cause a large upward shift in RMSE (e.g., compare Figures 18a and18c). This implies that our proposed inversion approach can have troubles to converge when applied to rich datasets with low error (e.g., ≤ 0.01 m for measured piezometric heads). We expect that a yet to be developed "localized" DR approach would perform better.

There are other aspects that would benefit from improvements. First, our DR approach is not totally accurate for direct conditioning. Second, several tens of thousands of MPS-generated training models are used to build it. This might not be practical when MPS-generation is computationally intensive. Third, 3D model generation capabilities could be enhanced by considering 3D convolutional layers instead of tricking 2D convolutional layers as done herein. Fourth, using a different reconstruction loss function than equation ( 8) the approach could be extended to categorical TIs with more than two facies and to continuous TIs. Lastly, it may be considered as a disadvantage that our approach produces continuous realizations and thus requires thresholding to create categorical fields.

Because of computational constraints, the impact of using a different number of dimensions for z on geostatistical simulation quality was not studied extensively. Yet limited testing revealed that for the considered 2D channelized TI and associated model domain size using a 25-dimensional z does not produce sufficient variability in the realizations (not shown). In contrast, compared to a 50-dimensional z using 100 dimensions induces fewer broken channels and thus provides slightly better results (not shown). For the considered 3D Maules Creek aquifer TI, using a 25-dimensional z was found to perform equally well as a 50-dimensional z (not shown). These findings warrant further investigations.

We would like to stress that if no hard thresholding is done at the end of our model generation process, then the derivatives of the pixels/voxels of the generated model, X, with respect to the elements of the associated low-dimensional code, z, have a (complicated) analytical solution that can be calculated by THEANO using auto-differentiation.

In other words, the sensitivity (Jacobian) matrix dX dz is immediately available. This may prove useful for gradient-based inversion using a forward solver equipped with an adjoint model.

In addition, it is worth noting that using a deep generative model to compress high dimensional quantities into a low-dimensional standard normal (or uniform) base is not only useful for inversion purposes but also holds promise for the so-called predictionfocused approach (e.g., [START_REF] Satija | Direct forecasting of reservoir performance using production data without history matching[END_REF]. The latter seeks to build a statistical model that links past data variables and prediction data variables, thereby enabling to make predictions without inverting for material properties. Lastly, note that our approach is not limited to the MPS framework but could potentially be applied to any prior with discrete structure.

Conclusion

This paper presents a novel parametric low-dimensional representation of complex binary geologic media that can be used for fast sampling of complex geostatistical models, which is necessary for efficient probabilistic inversion. A deep neural network of the variational autoencoder type is used to define a low-dimensional manifold that is restricted to model realizations that agree well with the training set. By comparing other parametric dimensionality reduction (DR) techniques for unconditional geostatistical simulation of a channelized prior model, we find that our representation is superior to the PCA, optimization-PCA (OPCA) and discrete cosine transform (DCT). Since constructing the parameterization requires a training set of several tens of thousands of model realizations obtained by multiple point statistics simulation, our DR approach is specifically designed for probabilistic (or deterministic) inversion rather than for unconditional (or point-conditioned) geostatistical simulation. Synthetic inversion case studies involving 2D steady flow and 3D transient hydraulic tomography are used to illustrate the effectiveness of the proposed DR-based probabilistic inversion approach. For the 2D case study, our new approach provides better results than current state-of-the-art inversion by sequential geostatistical resampling (SGR). Inversion results for the 3D application are also encouraging. Future work will focus on improvements such as alleviating the need of using a large training model set to build the parameterization, extension to multicategorical and continuous variables and defining a low-dimensional representation that is local rather than global. We expect that a local representation in which a given lowdimensional variable only influences a sub-region of the generated model realization will provide a much improved inversion performance.

  displays the corresponding CF metrics. It is observed that the DR-based realizations have similar CF curves as the training set. In terms of space of uncertainty, the ratio of the d JS (equation (12)) of the DR-based realizations to the d JS of the training images, d JS (DR) d JS (T R) , is 1.03. Thus, the DR-based realizations show slightly more between-realization variability than the training set. Figure 7 illustrates corresponding results obtained by the PCA, OPCA and DCT dimensionality reduction methods for the same training set. The PCA and OPCA realizations were derived following Vo and Durlofsky (2014) using 70 random variables.

Figure 5 :Figure 6 :

 56 Figure 5: Unconditional model realizations for the 2D binary channelized aquifer TI: (a -d) four (randomly chosen) model realizations from the DS-based training set and (e-l) eight (randomly chosen) model realizations generated with our DR approach that is based on a deep neural network.

Figure 7 :

 7 Figure 7: Unconditional model realizations for the 2D binary channelized aquifer TI: (ad) four (randomly chosen) PCA-based model realizations and (e-h) corresponding OPCA realizations together with (i-l) four (randomly chosen) DCT-based model realizations.

Figure 8 :Figure 9 :

 89 Figure 8: Conditional model realizations for the 2D binary channelized aquifer TI: (a -d) four (randomly chosen) conditioned model realizations from the DS-based training set and (e-l) eight (randomly chosen) conditioned model realizations generated with our DR approach that is based on a deep neural network.

Figure 10 :

 10 Figure 10: Unconditional model realizations for the 3D Maules Creek binary aquifer TI: (a -d) four (randomly chosen) model realizations from the MPS-based training set and (e-l) five (randomly chosen) model realizations generated with our DR approach that is based on a deep neural network.

Figure 11 :Figure 12 :

 1112 Figure 11: Cluster or connectivity function (CF) for the 3D case study involving the Maules Creek binary aquifer without direct conditioning data (Figure 10). The blue lines denote the values associated with the training set. The solid blue line indicates the mean while the 2 dashed lines represent the minimum and maximum values at each lag. The green solid lines represent the 100 realizations generated by our DR approach. The red solid line is the mean of these DR-based realizations. The CF is calculated for each facies along directions. The x, y and z symbols signify the x, y and z-axes. The d xy , d yz and d xz symbols represent the diagonal direction formed by the 45°angle from the x-axis in the xy plane, from the y-axis in the yz plane and from the x-axis in the xz plane, respectively.

Figure 13 :

 13 Figure 13: Trace plot of the sampled RMSE values by (a) the 8 Markov chains evolved by DREAM (ZS) (colored lines) for the first test (DR-generated true model) of inverse case study 1, (b) the 8 Markov chains evolved by DREAM (ZS) (colored lines) for the second test (DS-generated true model) of inverse case study 1 and (c) the 16 Markov chains evolved by DREAM (ZS) (colored lines) for the inverse case study 2. The dashed black line denotes the true RMSE value.

  Figure 14: (a) True model and (b-i) states of the 8 Markov chain evolved by DREAM (ZS) after 25,000 iterations per chain for the first test (DR-generated true model) of inverse case study 1. The red circle and blue crosses in subfigure (a) mark the location of the pumping well and the piezometers, respectively.

  Figure 15: (a) True model and (b-i) states of the 8 Markov chains evolved by DREAM (ZS) after 25,000 iterations per chain for the second test (MPS-generated true model) of inverse case study 1. The red circle and blue crosses in subfigure (a) mark the location of the pumping well and the piezometers, respectively.

  Figure 16: (a) True model and (b-i) states of the 8 independent SGR trials after 25,000 iterations per trial for the second test (MPS-generated true model) of inverse case study 1. The red circle and blue crosses in subfigure (a) mark the location of the pumping well and the piezometers, respectively.

Figure 17

 17 Figure 17: a) True model and (b-i) states of the first 8 (out of 16) Markov chains evolved by DREAM (ZS) after 60,000 iterations per chain for the inverse case study 2.

Figure 18 :

 18 Figure 18: Model realizations and associated RMSE obtained when modifying a single low-dimensional dimension, DR 1 , of a posterior solution for case study 1: (a) reference (posterior) realization, (b) realization and RMSE corresponding to the multiplication of DR 1 by 0.5, (c) realization and RMSE corresponding to the multiplication of DR 1 by 0.1 and (d) realization and RMSE corresponding to the multiplication of DR 1 by -1.

For instance:http://deeplearning.net/tutorial/ and http://cs231n.github.io/ convolutional-networks/.

Acknowledgments

Python codes of the proposed dimensionality reduction (DR) and DR-based inversion approaches are available from the first author (and can be downloaded from https: //github.com/elaloy). This work was partially supported by Agence Nationale de la Recherche through the grant ANR-16-CE23-0006 Deep in France. The third author is a Research Associate with the Belgian F.R.S.-F.N.R.S. We thank Mike Swarbrick Jones for sharing his codes (https://github.com/mikesj-public), one of which served as a starting point for our devised VAECNN. We are grateful to Hai Xuan Vo and Louis Durlofsky for sharing their OPCA code, and to Xue Li and her co-authors for sharing their GC MPS algorithm. A temporary academic license of the DeeSse MPS code can be obtained upon request to one of its developers (Grégoire Mariethoz, Philippe Renard, Julien Straubhaar). We thank Laurent Lemmens for sharing his MPH and CF calculation routines. Finally, we would like to thank the two anonymous referees for their useful comments.