DeepFlow
History Matching in the Space of Deep Generative Models
Abstract
The calibration of a reservoir model with observed transient data of fluid pressures and rates is a key task in obtaining a predictive model of the flow and transport behaviour of the earth’s subsurface. The model calibration task, commonly referred to as "history matching", can be formalised as an illposed inverse problem where we aim to find the underlying spatial distribution of petrophysical properties that explain the observed dynamic data. We use a generative adversarial network pretrained on geostatistical objectbased models to represent the distribution of rock properties for a synthetic model of a hydrocarbon reservoir. The dynamic behaviour of the reservoir fluids is modelled using a transient twophase incompressible Darcy formulation. We invert for the underlying reservoir properties by first modeling property distributions using the pretrained generative model then using the adjoint equations of the forward problem to perform gradient descent on the latent variables that control the output of the generative model. In addition to the dynamic observation data, we include well rocktype constraints by introducing an additional objective function. Our contribution shows that for a synthetic test case, we are able to obtain solutions to the inverse problem by optimising in the latent variable space of a deep generative model, given a set of transient observations of a nonlinear forward problem.
DeepFlow
History Matching in the Space of Deep Generative Models
Lukas Mosser Department of Earth Science and Engineering Imperial College London lukas.mosser15@imperial.ac.uk Olivier Dubrule Department of Earth Science and Engineering Imperial College London o.dubrule@imperial.ac.uk Martin J. Blunt Department of Earth Science and Engineering Imperial College London m.blunt@imperial.ac.uk
noticebox[b]Preprint. Work in progress.\end@float
1 Introduction
Understanding the flow of fluids within the earth’s subsurface is a key issue in many practical applications such as the understanding of how pollutants in aqueous phases are transported and affect the environment (Yong et al., 1992; Bear & Cheng, 2010), hydrogeological considerations (Domenico et al., 1998; Fetter, 2018), the formation of mineral deposits (Garven, 1985; Sibson et al., 1975), carboncapture and sequestration (CCS) (Lackner, 2003; Holloway, 2005), and the migration and production of hydrocarbons from subsurface reservoirs (England et al., 1987). When modeling subsurface architectures a spatial gridded representation is usually built for modeling physical flow processes (Aziz, 1993). The flow of fluids, including aqueous and nonaqueous fluids such as supercritical is simulated by solving a set of partial differential equations that define the transient behaviour of the flowing liquid phases (Gerritsen & Durlofsky, 2005).
In many cases the ability for fluids to flow within a porous material is dominated by two key properties: the effective porosity defines how much pore space is available to flow as a proportion of the bulk volume of the material (Bear, 2013). The permeability is a property that is proportional to the flow conductance (Darcy, 1856; Muskat, 1938). We refer to Blunt et al. (2013) for a detailed review of the principles of multiphasefluid flow within porous media. Obtaining direct measurements of these two properties, permeability and porosity, within the earth’s structures is difficult and in many cases these measurements are sparse and obtained at various spatial scales, making their integration within the modeling workflow difficult (Wu et al., 2002; Farmer, 2002).
Nevertheless, representing the spatial distribution of the properties of porous media within the earth’s subsurface is necessary to build reliable predictive models. In applications such as the injection of supercritical into hosting rock formations for longterm storage (Chadwick & Noy, 2010; Martens et al., 2012) or the production of hydrocarbons and water from subsurface reservoirs we obtain, in addition to well and seismic data, measurements of the variables dependent on these flow properties ( and ) such as fluid phase pressures, produced and injected volumes over time (Muskat, 1949). Calibration of a spatial model of the petrophysical and structural features of the earth’s subsurface with available physical measurements is commonly referred to as "history matching" (Oliver & Chen, 2011), and more generally as data assimilation (Reichle et al., 2002). Determining the underlying spatial distributed parameters from observed parameters such as fluid pressures and rates, can be formalised as an an illposed inverse problem.
For an inverse problem to be wellposed it must satisfy three criteria (Tikhonov & Arsenin, 1977)

Existence: for every set of observed parameters there must exist a solution.

Uniqueness: the obtained solution has to be unique.

Stability: the solution is stable with regards to small changes in the observed parameters i.e. a small change in the observed parameters induces only a small change in the solution parameters
if these three criteria are not fulfilled the problem is considered illposed. It has been shown that inverse problems with spatially distributed parameters common in groundwater and reservoir flow problems are often illposed because they don’t fulfill the uniqueness or stability criterion (Kravaris & Seinfeld, 1985). Chavent (1979) provides a detailed analysis of the existence criterion for distributed parameter systems of the heat, transport, and wave equations.
Reducing the number of parameters used to represent the solution space allows the inverse problem to be welldefined. Coats et al. (1970) assumed a constant value for large spatial regions of the underlying parameters (Cooley & Sinclair, 1976; Yeh & Tauxe, 1971). Leo et al. (1986) used Bsplines to represent spatial distributions of porosity and permeability. Recently, a number of approaches for model parameterisation have been developed based on principal component analysis (PCA) (Vo & Durlofsky, 2014) and improvements thereof using kernel and optimisationbased PCA. Demyanov et al. (2011) use kernel learning to find a parameterisation of the petrophysical parameter space.
In a Bayesian setting of the inverse problem (Tarantola, 2005) the gradual deformation method applies a continuous perturbation of a set of model parameters obtained from a prior distribution to honor the observed data of the forward problem (Roggero & Hu, 1998; Caers, 2007). Bayesian approaches to solve the reservoir history matching problems, such as the Ensemble Kalman Filter (EnKF) (Evensen, 1994, 2003) have been applied in numerous historymatching case studies (Lorentzen et al., 2005). Emerick & Reynolds (2012) present an ensemble method based on EnKF, that shows improved performance in the case of highly nonlinear problems.
Recent developments in deep learning (LeCun et al., 2015; Goodfellow et al., 2016) have motivated the use of neuralnetworkbased parameter representations. Canchumuni et al. (2018) parameterised geological training images using a Deep Belief Network (DBN) (Hinton et al., 2006; Hinton & Salakhutdinov, 2006) and used an ensemble smoother with multiple data assimilation (ESMDA) (Emerick & Reynolds, 2013) to perform history matching. Variational autoencoders (Laloy et al., 2017) and spatial generative adversarial networks (SGAN) (Jetchev et al., 2016; Laloy et al., 2018) have been used to generate a set of channelised training images to perform probabilistic inversion using adaptive MarkovChain MonteCarlo methods (Laloy & Vrugt, 2012). Chan & Elsheikh (2018a) used generative adversarial networks trained using patchbased kernel discrepancy measures to synthesise geological models from training images. Realisations conditional to spatial observations of the underlying model parameters can be obtained by gradientbased image inpainting methods (Yeh et al., 2016; Mosser et al., 2018a; Dupont et al., 2018) or direct learning of conditional image distributions using autoregressive generative models such as PixelCNNs (den Oord et al., 2016; Dupont & Suresha, 2018).
Adjointstate methods allow gradients of the forward problem with respect to the solution parameters to be obtained. These gradients can then be further backpropagated through a deep generative models such as generative adversarial networks (Goodfellow et al., 2014; Kingma & Welling, 2013; Richardson, 2018; Mosser et al., 2018b) to perform inversion in the space of latent variable models and have been used for probabilistic inversion in the case of acoustic waveform inversion using a Metropolisadjusted Langevin algorithm (MALA) (Roberts & Stramer, 2002; Mosser et al., 2018b). Numerous studies have evaluated the use of deep generative models such as GANs for linear inverse problems (Shah & Hegde, 2018) in the context of compressed sensing (Bora et al., 2018; Mardani et al., 2019) and have analysed the convergence of inversion algorithms that use GANs (Bora et al., 2017; Shah & Hegde, 2018).
A number of recent approaches have combined methods developed in the context of deep neural networks with traditional approaches to solve (illposed) inverse problems such as learning the regularisation itself (Li et al., 2018; Lunz et al., 2018), using iterative deep neural networks (Adler & Öktem, 2017), or introducing deep neural networks that account for nonlinear forward operators in reconstruction algorithms (Adler & Öktem, 2018).
Following the work of Mosser et al. (2018b) we parameterise a family of geological models and associated permeability and porosity values using a GAN. We solve the transient twophase immiscible Darcy flow equations without gravity and under isothermal conditions using a finitedifference approach and obtain gradients of the mismatch between observed and simulated data with respect to the gridblock permeability and porosity using the adjointstate method.
Starting from random initial locations within the underlying Gaussian latent space controlling the GAN’s output, we use gradientdescent to minimise the misfit between the flow data and the forward model output by modifying the set of latent variables. A regularisation scheme of the optimisation problem is derived from a Bayesian inversion setting. Based on a twodimensional synthetic case of a channelised reservoir system we show that gradientbased inversion using GANs as a model parameterisation can be achieved, honoring the observed physical quantities of the nonlinear forward problem as well as borehole measurements of the underlying petrophysical quantities. The trained models and code are available as opensource software^{1}^{1}1https://github.com/LukasMosser/DeepFlow.
2 Methodology
We consider the problem of determining the flow behaviour of a subsurface hydrocarbon reservoir from which oil and water are being produced from a well. To displace the oil phase we inject water at a constant injection rate in another well and produce oil by maintaining a constant bottomhole pressure. The forward problem is therefore defined as solving the twophase (oil and water) Darcy flow equations with a slightly compressible oil phase with no gravity effects and under isothermal conditions. The observed variables of the system are the injection pressure at the location of the injection well as a function of time , as well as the produced volumes of water and oil at the production well .
(1) 
We consider all fluid parameters of the forward problem, the oil and water viscosity , density , compressibility , as well as the relativepermeability behaviour, represented by a BrooksCorey model, as constants (Blunt, 2017). A table of fluid properties is provided in the appendix (see table 2). Therefore, the set of model parameters are the spatial distribution of the reservoir rock type indicator function , the singlephase permeability , and porosity .
(2) 
We assume a binary set of reservoir rocktypes where corresponds to a low permeability shale and indicates the high permeability sandstone of a riverchannel.
Given the nonlinear forward operator this allows us to define the forward problem as
(3a)  
(3b) 
We represent the distribution of possible model parameters by a (deep) generative latentvariable model
(4) 
with latent variables sampled from a Gaussian distribution (Eq. 4). Here we assume that the generative model has been trained prior to the inversion process, and is therefore constant.
Furthermore we can write:
(5) 
and from Bayes theorem
(6) 
where for the observed flow data the likelihood function is a Gaussian with mean and variance . We have assumed a Bernoulli likelihood for the observed facies at the wells (Eq. 2).
Furthermore, we assume conditional independence of the likelihood of the observed flow and facies data at the wells
(7) 
In the case of GANs only individual samples can be obtained and therefore estimation of the posterior is restricted to pointestimates. We derive a regularisation scheme for the history matching problem from finding point estimates of the maxima of the posterior distribution (MAP)
(8a)  
(8b)  
(8c)  
(8d)  
(8e)  
where the squarederror norm of the difference between the observed data and the modelled data is a result of the assumption of a Gaussian likelihood in equation 3 and the well loss corresponds to the binary crossentropy between the observed facies indicator and the facies probability of each gridblock at the wells.
The derived regularisation scheme can be framed as a minimisation problem over the latent variables i.e. we minimise Eq. 8 by modifying the latent variables . Specifically we use a squarederror functional that measures the difference between the observed pressures and rates at the well locations averaged over the total duration of the observed data
(9) 
where and are the locations of the production and injection wells respectively and measurement uncertainties have been quantified by the standard deviation of the oil and water rate
(10) 
and the standard deviation of the pressure measurements
(11) 
used to generate the set of observed pressure and rate data.
Using the adjointstate method (Plessix, 2006; Suwartadi et al., 2012; Krogstad et al., 2015) and traditional neural network backpropagation (Rumelhart et al., 1988) we can obtain gradients of the error functional (Eq. 9) with respect to the latent variables.
We therefore seek to minimise the total loss
(12) 
consisting of the flow loss (Eq. 9), the welldata loss , and the prior loss (Eq. 8). We note that this approach is similar to the approach by Creswell & Bharath (2018) which includes an additional loss term corresponding to the prior distribution of the latent variables when inverting the generator function of a GAN. The system of partial differential equations of the nonlinear forward problem as well as the likelihood function (Eq. 9) have been previously implemented in the opensource MATLAB/GNUOctave reservoir simulation framework MRST (Lie, 2019). An overview of the endtoend coupling between the latent variables, the deep generative model , and the numerical forward model is shown in figure 1.
To evaluate the proposed approach we create synthetic twodimensional vertical crosssections of a stacked riverchannel system using an objectbased approach. Channelbodies are represented by halfcircles with centers following uniform distributions in space. The full objectbased modeling workflow is detailed in Mosser et al. (2018b). For the case of reservoir history matching (Eq. 3) each objectbased model is associated with the three model parameters following equation 2 (Fig. 2).
We train a deep generative adversarial network (GAN) (Goodfellow et al., 2014) on a set of objectbased realisations to learn the generative model (Eq. 4) prior to obtaining MAP samples (Eq. 8) by optimising equation 12. Specifically, we train a socalled WassersteinGAN (Arjovsky et al., 2017; Gulrajani et al., 2017), that minimises the Wasserstein distance between generated and real probability distributions. To stabilise training of the WassersteinGAN we use the Lipschitz penalty introduced by Petzka et al. (2017). A detailed overview of the training procedure for the deepgenerative model is presented in Mosser et al. (2018b). Network architectures for the generator and discriminator are detailed in appendix table 1.
The generator network is a deep convolutional neural network that maps a 100dimensional noisesample from a standardised multiGaussian distribution to a spatial distribution of the model parameters (Eq. 2) with 128 pixels in the x and 64 pixels in the zdirection
(13) 
where the three output features correspond to the gridblock rocktype probability, permeability and porosity (Eq. 2). Due to the continuous nature of the generative networks sigmoid activation function, the output of the generative network corresponding to a representation of the reservoir rocktype can be interpreted as a rocktype probability. The gridblock permeability and porosity values are obtained by transforming the generated rocktype probability from the interval to a permeability by a linear transformation (see appendix table 1(a)) to rescale the output of the generative model to a range of values that correspond to admissible petrophysical values for sandstone and shale. The finitedifference approximation of the forward model is computed on the regular grid given by the pixelbased output of the generative model and hence each pixelbased property corresponds to the gridblock property used in the numerical evaluation of the forward problem.
We have selected a reference realisation from a testset of objectbased realisations that were generated independently of the models used to train the generative network . The model shown in figure 2 was used to generate a set of synthetic observed data that served as a reference case used to validate the proposed inversion approach. Gaussian noise was added to the production rates and pressures based on assumed measurement uncertainties (Eq. 1011). The transient pressure and rates of oil and water were recorded for a duration of 600 days and the full history was used in the gradientbased optimisation outlined in equation 8. To minimise the total loss (Eq. 12) we used the ADAM optimiser (Kingma & Ba, 2014) using a fixed step size and parameters (, , ). The simulation parameters used for the numerical evaluation of the forward problem are provided in appendix table 2.
3 Results
Four history matching scenarios were considered:

Scenario 1: Samples of the unconditional prior were obtained () and a single evaluation of the forward problem was performed. This allowed us to evaluate the prior distribution of the dynamic flow behaviour of the underlying generative model (Fig. 3a).

Scenario 2: Inversion was performed by considering only the well data loss at the injection and production wells and the prior loss to honor the prior distribution of the latent variables (Fig. 3b).

Scenario 3: We performed inversion by considering only the flowbased loss due to the mismatch between simulated and observed dynamic twophase flow data, as well as the prior loss on the latent variables (Fig. 3c).
For scenario 2 to 4, we performed inversions, starting from the same set of initial random samples of latent variables . Each optimisation procedure was run for a total of 500 ADAM optimisation steps. For scenarios 3 and 4, we choose the iteration with the lowest total loss as the output of each run. Early stopping was performed in scenario 2 where only the well data loss was considered when the facies accuracy at the well reached 100.
Figure 3 shows a comparison of these four scenarios in terms of their dynamic flow behaviour i.e. pressures and rates at production and injection wells as a function of the simulation time, as well as the mean and standard deviation of the thresholded rocktype probability maps. The threshold for rocktype probability maps was set at 0.5 to distinguish between sand (1) and shale (0).
The prior distribution of the observed variables (Scenario 1, Fig. 3a, right) shows a large spread in the distribution of the observed oilwater production rates and injection pressure. The observed liquid rates of the reference case (red) lie at the edge of the unconditional data distribution. The pressure behaviour of the prior models varies across 2 orders of magnitude in relation with the pressure behaviour of the reference case at the low end of the pressure versus time distribution.
The unconditional prior samples of the model parameters, shown here as the rocktype indicator function , are shown as crosssections of the mean and standard deviation of the models (Scenario 1, Fig. 3a, left). As shown by the standard deviation map, nearly all values are close to a value of 0.5 as expected from an indicator function of mean 0.5.
Considering the likelihood of the gridblock rocktype indicator function at the wells only (Scenario 2) where the binary crossentropy loss at the injection and production well , as well as the prior loss were minimised, we find a large spread in the production behaviour of the generated samples closely resembling the unconditional prior distribution (Scenario 2, Fig. 3b, right). Inspecting the mean and standard deviation maps of the rocktype indicator function of the inverted samples matching the well data only, we find that these samples match the well data while their flow behaviour is consistent with the unconditional prior. This means that constraining only by well data has very little effect on the facies model connectivity between the wells. A lateral effect of the gridblock well data can be observed comparing the outline of the mean rocktype maps and the red contours of the ground truth model (Fig. 2) used to generate the reference observed data and .
When only the flow loss was considered in the optimisation, the production data of the generated realisations were tightly constrained to the observed liquid rates and pressures (Scenario 3, Fig. 3c, right). The mean and standard deviation maps for the obtained distribution of the rocktype indicator functions do not indicate any clear structural features and resemble the prior maps shown in the first row of figure 3. We have obtained samples of model parameters that all honor the observed dynamic flow data , but which are not anchored at the wells by the observed well rock types.
Combining the flow loss , the well data loss and prior loss as the total loss (Eq. 12) (Scenario 4, Fig. 3d) we obtain samples that match the observed dynamic behaviour closer than scenario three, which only considered the flow loss. In addition to honoring the dynamic observed data , the facies data of the obtained samples at the wells match the observed gridblock data at the production and injection well . This is indicated by the low standard deviation map of the inverted models (Scenario 4, Fig. 3d, left).
Figure 4 compares the distribution of the total loss (, Eq. 12) for the models obtained from the unconditional prior distribution i.e. sampling a set of model parameters (Eq. 2) from the generative model and solving the forward problem once, (Eq. 3) and the models of scenario 3 (flow and prior loss only) and 4 (flow, well and prior loss). For both considered scenarios (Fig. 4, a and c respectively) the distribution of the total loss is shifted to lower values when compared to the total loss distribution of the models obtained by sampling from the prior distribution. The distribution of well rock type accuracy shows that history matching only flow data (Fig. 4 b) does not improve the rock type accuracy. When the welldata loss is optimised in addition to the flowloss (Fig. 4 d), we observe that the obtained models closely match the rock type at the wells in addition to the dynamic flow data.
Due to the high computational cost of the forward model, reaching a good match with the observed data, within a reasonable number of forwardmodel evaluations, is necessary. Figure 9 shows the distribution of the number of optimisation iterations of equation 8 necessary to reach a given threshold value of the total loss for the case where all three objective functions are used in the history matching process (Scenario 4). All inverted models achieve a totalloss less than (Fig. 9 c) which corresponds to the low end of the loss values obtained by sampling from the prior distribution only (Fig. 4 c).
In figure 3 we have shown that including the flow loss leads to high variance in the generated distribution of the model properties while honoring the observed flow data (Scenario 3, Fig. 3c, right), while only constraining samples to honor the well data reduces the variance near the wells, but does not constrain to observed dynamic data (Scenario 2, Fig. 3b, right). A smaller injection pressure observed for the samples obtained from scenario 3 compared to those that match only the wells in scenario 2 would indicate a better connected system of channel bodies. After closer inspection of the inverted samples it was found that nearly all models obtained from scenario 3 and 4 are connected (Fig. 5 and appendix Fig. 1314). We evaluate for all four scenarios whether any of the clusters of riverchannel bodies connects the injection and production wells. Less than half of the models for all scenarios showed a connection between injector and producer wells (Appendix Fig. 5  black). After Boolean dilation of the largest cluster and performing the connectivity analysis again, half of the samples that honor the well data only (Scenario 2) show connectivity, whereas all models obtained that were optimised to honor flow (Scenario 3), as well as flow and well data (Scenario 4) form a connected cluster between the injector and producer well (Appendix Fig. 5  red). A possible explanation of this result is that due to the finitedifference approximation (twopoint flux approximation) used to solve the forwardproblem the transmissibilities between two neighboring grid blocks are computed as a weighted average of the permeability values i.e. a single grid block spacing between two highpermeable channelbodies does not provide an effective flow barrier.
A single model obtained from optimisation using flow and well losses, while honoring the prior distribution of the latentvariables (Scenario 4), was chosen as an example to highlight the optimisation process for a single realisation. This specific example was chosen as it showed the largest reduction in the total loss from the initial starting model to the iteration with the lowest total loss out of 500 optimisation steps. Figure 6 shows the evolution of the total loss (black), the accuracy in matching the rocktype indicator function at the well locations (red), and the prior loss evolution. A highly nonlinear optimisation with many local minima can be observed. The accuracy of matching the well data can be observed to quickly reach values above 90% and staying at 100% for the entirety of the optimisation process. The smallest total loss for this model is achieved after 214 iterations.
For the chosen example model, we have recorded the observed dynamic flow behaviour as well as (left column) the distribution of the flow parameters at intermediate steps in the optimisation procedure, where the total loss achieved a new local minimum. Figure 7 shows the evolution of the dynamic flow behaviour in the course of performing the inversion. The model initially shows a very high injection pressure (Fig. 7 c blue), which is significantly and continuously reduced to values close to the observed pressure data (red). Similar behaviour is observed for the liquid water and oil rates at the production wells (Fig. 7 ab). The solution found by the optimisation strategy (Eq. 8) with the lowest error value (Fig. 7 ac, black) closely matches the observed dynamic production data and the rocktype indicator function at the wells.
The intermediate spatial distributions of the rocktype indicator for each of the obtained local minima in the optimisation process are shown in figure 8. Initially the channelbodies were completely disconnected (Fig. 8 a). As the optimisation process progresses more channelbodies were introduced gradually (Fig. 8 be) and form a spanning cluster that matches the gridblock data at the wells (Fig. 8 f).
4 Discussion
We have shown that solutions of the illposed inverse problem of reservoir history matching can be obtained by inversion in the space of a deep generative adversarial network for a synthetic case. The obtained model parameter distributions match the dynamic observed flow data as well as honor the gridblock properties at the wells with adequate solutions found on average in less than 100 optimisation iterations. To evaluate whether the obtained samples of the posterior distribution are local maxima we followed the approach by Zhang et al. (2019) and performed a cyclical spherical linear interpolation (White, 2016) in the latent space between three obtained realisations. For each interpolated latent vector we have solved the forward model and evaluated the objective functions. Figures 8(a) and 8(b) indicate that the posterior distribution is highly multimodal, where each of the obtained realisations lie at a minimum of the total loss i.e. near a local maximum of the posterior.
While the obtained samples match the flow and well data it is important to note that in a highdimensional setting, finding the maxima of the posterior may not be a sufficient description of the full posterior distribution and this needs to be considered when using the obtained historymatched ensemble of reservoir models for forecasts of future reservoir performance. To estimate posterior parameters such as the cumulative production, a harmonic mean averaging approach (Green, 1995; Raftery et al., 2006) could be used where individual samples are weighted by their corresponding likelihood (Zhang et al., 2019). Furthermore, a number of open challenges are still to be addressed in the context of solving inverse problems with deep generative models such as GANs.
The presented synthetic case is a simplified representation of a system of riverchannel bodies. The features encountered in real reservoir systems are more diverse and heterogeneous, and it is possible that many other reservoir features should be included in the training set of generative models. The types of models that may be needed include threedimensional multirocktype models i.e. more than two rocktypes in a generated reservoir model. Furthermore, it may be of interest to combine a number of datasets into a single generative model that is conditioned on the type of reservoir architecture , where is a vector of conditional variables that influence the output of the deep generative model using discrete conditionals such as depositional environment e.g. fluvial versus lacustrine environments or using continuous conditioning variables such as reservoir nettogross. This type of conditional generative model could be facilitated by making use of architectures such as that proposed by Miyato & Koyama (2018). Empirically, it has been shown that introducing classconditional variables to GANs, socalled cGANs, result in improved training and distributional representation behaviour compared to their unconditional counterparts (Gauthier, 2014). This may speak for training a single general reservoir model generator compared to creating a library of standalone pretrained generator networks for each depositional environment.
An important aspect to be considered here is the evaluation of the generated models in terms of model quality, modecollapse, and representation of spatial statistics as highlighted in Mosser et al. (2017) and Laloy et al. (2018), where twopoint correlation functions are used to quality control the individual generated reservoir models. The use of established GAN quality control measures such as Inception Score (IS) (Salimans et al., 2016) or Frechet Inception Distance (FID) (Heusel et al., 2017) may be challenging due to the induced distribution shift between ImageNet (Deng et al., 2009) pretrained networks and the features observed in typical reservoir architectures. Kernelbased methods such as the maximum mean discrepancy (MMD) have already been successfully applied to train GANs on reservoir models (Chan & Elsheikh, 2018b) which could be repurposed to quality control generated reservoir models.
While the synthetic case presented herein has shown that a gradientbased approach to solving the illposed inverse problem is feasible, other strategies such as using MarkovChain MonteCarlo methods as demonstrated by Laloy et al. (2017), may lead to successful inversion results, but at high computational cost. Combinations with other established inversion techniques, such as Ensemble KalmanFilters (Gu & Oliver, 2005; Aanonsen et al., 2009) or Ensemble Smoothers (ESMDA) (Emerick & Reynolds, 2013) should also be evaluated. The presented approach, together with a set of suitable generative models should be applied to benchmark studies such as the PUNQ or SPE comparative solutions datasets to investigate the applicability of inversion with deep generative models in the presence of measurement errors and noise. While the present study only considers inverting for one model at a time i.e. using a batchsize of one, it is possible to perform inversion using a batchsize greater than one, as in Creswell & Bharath (2018) who show that the gradient contributions to the latent variables can be computed independently for each element of a batch of latentvariables and hence for a set of model parameters . By computing the numerical solution forward problem and the associated adjoint for each set of modelparameters within a batch, the proposed method could be integrated into existing workflows for ensemble reservoir modeling.
Using deep latent variable models, such as GANs, to represent the space of solutions for illposed inverse problems opens up interesting avenues with regards to the theoretical basis of the illposed nature of these problems. Investigation of the existence and stability criteria (Sec. 1) (Tikhonov & Arsenin, 1977) of deep latent variable inversion schemes could lead to new promising theoretical insights into the convergence behaviour of inversion methods using deep generative models. Furthermore, the geometry of the loss landscape and representation manifold within the space of latent variables may allow for faster convergence and higher quality inversion results. Chen et al. (2017) and Arvanitidis et al. (2017) have investigated the geometric structure of the manifold embedded in the latentspace of deep generative models and Shao et al. (2018) have shown that the manifolds obtained from VAEs are nonlinear with nearzero curvature and and hence linear paths in latent space correspond to geodesics on the embedded manifold. This indicates that considering the geometry of the learned manifolds could be accounted for in operations performed in the latent space of deep generative models.
We have shown that GANs may be used as generative models for parameterisation of geological models, and can provide a solution space for illposed inverse problems. Nevertheless, GANs due to their challenging training, evaluation and complex latenttoparameterspace relationship may not be the optimal reparameterisation choice in the case of illposed inverse problems. Deep generative models such as variational (Kingma & Welling, 2013) and disentangling autoencoders (Higgins et al., 2017; Burgess et al., 2018), as well as flowbased generative models, such as RealNVP (Dinh et al., 2016) or GLOW (Kingma & Dhariwal, 2018), that provide a bijective and invertible mapping may be better suited for the inversion task. This has been demonstrated by Ardizzone et al. (2018) for illposed inverse problems.
5 Conclusions
We have presented an application of deep generative models in the context of the illposed history matching inverse problem. Based on a twodimensional synthetic crosssection of a riverchannel system we have trained a GAN to represent the prior distribution of the subsurface properties , permeability , porosity and the rocktype indicator . We find solutions to the twophase slightlycompressible Darcy flow problem of oil and water systems commonly used to describe the transient flow behaviour in hydrocarbon reservoirs. By incorporating a finitedifferencebased numerical simulator with adjointstate capabilities (Krogstad et al., 2015) in an endtoend differentiable framework, we perform inversion using gradientbased optimisation of the mismatch between observed dynamic data (Eq. 9) and gridblockscale well data (Fig. 3), while honoring the prior distribution of the latent variables (Eq. 12). By using a momentumaccelerated firstorder gradient descent scheme (Kingma & Ba, 2014), the method converges to a solution despite the highly nonlinear and nonconvex loss landscape (Fig. 6). Future work will focus on applications to real reservoirs and history matching benchmark studies such as the PUNQ or SPE comparative solutions models, as well as evaluation of other deep latentvariable generative models and their theoretical benefits for finding solutions to illposed inverse problems.
Acknowledgments
O. Dubrule would like to thank Total S.A. for seconding him as visiting professor at Imperial College London.
References
 (1)
 Aanonsen et al. (2009) Aanonsen, S. I., Nævdal, G., Oliver, D. S., Reynolds, A. C. & Vallès, B. (2009), ‘The ensemble Kalman filter in reservoir engineering–a review’, SPE Journal 14(03), 393–412.
 Adler & Öktem (2017) Adler, J. & Öktem, O. (2017), ‘Solving illposed inverse problems using iterative deep neural networks’, Inverse Problems 33(12), 124007.
 Adler & Öktem (2018) Adler, J. & Öktem, O. (2018), ‘Learned primaldual reconstruction’, IEEE Transactions on Medical Imaging 37(6), 1322–1332.
 Ardizzone et al. (2018) Ardizzone, L., Kruse, J., Wirkert, S., Rahner, D., Pellegrini, E. W., Klessen, R. S., MaierHein, L., Rother, C. & Köthe, U. (2018), ‘Analyzing inverse problems with invertible neural networks’, arXiv preprint arXiv:1808.04730 .
 Arjovsky et al. (2017) Arjovsky, M., Chintala, S. & Bottou, L. (2017), ‘Wasserstein gan’, arXiv preprint arXiv:1701.07875 .
 Arvanitidis et al. (2017) Arvanitidis, G., Hansen, L. K. & Hauberg, S. (2017), ‘Latent space oddity: on the curvature of deep generative models’, arXiv preprint arXiv:1710.11379 .
 Aziz (1993) Aziz, K. (1993), ‘Reservoir simulation grids: opportunities and problems’, Journal of Petroleum Technology 45(07), 658–663.
 Bear (2013) Bear, J. (2013), Dynamics of fluids in porous media, Courier Corporation.
 Bear & Cheng (2010) Bear, J. & Cheng, A. H.D. (2010), Modeling groundwater flow and contaminant transport, Vol. 23, Springer Science & Business Media.
 Blunt (2017) Blunt, M. J. (2017), Multiphase flow in permeable media: A porescale perspective, Cambridge University Press.
 Blunt et al. (2013) Blunt, M. J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P., Paluszny, A. & Pentland, C. (2013), ‘Porescale imaging and modelling’, Advances in Water Resources 51, 197–216.
 Bora et al. (2017) Bora, A., Jalal, A., Price, E. & Dimakis, A. G. (2017), Compressed sensing using generative models, in ‘Proceedings of the 34th International Conference on Machine LearningVolume 70’, JMLR. org, pp. 537–546.
 Bora et al. (2018) Bora, A., Price, E. & Dimakis, A. G. (2018), AmbientGAN: Generative models from lossy measurements, in ‘International Conference on Learning Representations (ICLR)’.
 Burgess et al. (2018) Burgess, C. P., Higgins, I., Pal, A., Matthey, L., Watters, N., Desjardins, G. & Lerchner, A. (2018), ‘Understanding disentangling in betaVAE’, arXiv preprint arXiv:1804.03599 .
 Caers (2007) Caers, J. (2007), ‘Comparing the gradual deformation with the probability perturbation method for solving inverse problems’, Mathematical Geology 39(1), 27–52.
 Canchumuni et al. (2018) Canchumuni, S. W., Emerick, A. A. & Pacheco, M. A. C. (2018), History matching channelized facies models using ensemble smoother with a deep learning parameterization, in ‘ECMOR XVI16th European Conference on the Mathematics of Oil Recovery’.
 Chadwick & Noy (2010) Chadwick, R. A. & Noy, D. J. (2010), Historymatching flow simulations and timelapse seismic data from the Sleipner CO2 plume, in ‘Geological Society, London, Petroleum Geology Conference series’, Vol. 7, Geological Society of London, pp. 1171–1182.
 Chan & Elsheikh (2018a) Chan, S. & Elsheikh, A. H. (2018a), ‘Exemplarbased synthesis of geology using kernel discrepancies and generative neural networks’, arXiv preprint arXiv:1809.07748 .
 Chan & Elsheikh (2018b) Chan, S. & Elsheikh, A. H. (2018b), ‘Parametric generation of conditional geological realizations using generative neural networks’, arXiv preprint arXiv:1807.05207 .
 Chavent (1979) Chavent, G. (1979), ‘Identification of distributed parameter systems: about the output least square method, its implementation, and identifiability’, IFAC Proceedings Volumes 12(8), 85–97.
 Chen et al. (2017) Chen, N., Klushyn, A., Kurle, R., Jiang, X., Bayer, J. & van der Smagt, P. (2017), ‘Metrics for deep generative models’, arXiv preprint arXiv:1711.01204 .
 Coats et al. (1970) Coats, K. H., Dempsey, J. R. & Henderson, J. H. (1970), ‘A new technique for determining reservoir description from field performance data’, Society of Petroleum Engineers Journal 10(01), 66–74.
 Cooley & Sinclair (1976) Cooley, R. L. & Sinclair, P. J. (1976), ‘Uniqueness of a model of steadystate groundwater flow’, Journal of Hydrology 31(34), 245–269.
 Creswell & Bharath (2018) Creswell, A. & Bharath, A. A. (2018), ‘Inverting the generator of a generative adversarial network’, IEEE transactions on neural networks and learning systems .
 Darcy (1856) Darcy, H. P. G. (1856), Les Fontaines publiques de la ville de Dijon. Exposition et application des principes à suivre et des formules à employer dans les questions de distribution d’eau, etc, V. Dalamont.
 Demyanov et al. (2011) Demyanov, V., Foresti, L., Christie, M. A. & Kanevski, M. (2011), Reservoir modelling with feature selection: kernel learning approach, in ‘SPE Reservoir Simulation Symposium, SPE 141510MS’, Society of Petroleum Engineers.
 den Oord et al. (2016) den Oord, A., Kalchbrenner, N., Espeholt, L., Vinyals, O. & Graves, A. (2016), Conditional image generation with PixelCNN decoders, in ‘Advances in Neural Information Processing Systems’, pp. 4790–4798.
 Deng et al. (2009) Deng, J., Dong, W., Socher, R., Li, L.J., Li, K. & FeiFei, L. (2009), ‘Imagenet: A largescale hierarchical image database’.
 Dinh et al. (2016) Dinh, L., SohlDickstein, J. & Bengio, S. (2016), ‘Density estimation using real NVP’, arXiv preprint arXiv:1605.08803 .
 Domenico et al. (1998) Domenico, P. A., Schwartz, F. W. & Others (1998), Physical and chemical hydrogeology, Vol. 506, Wiley New York.
 Dupont & Suresha (2018) Dupont, E. & Suresha, S. (2018), ‘Probabilistic semantic inpainting with pixel constrained CNNs’, arXiv preprint arXiv:1810.03728 .
 Dupont et al. (2018) Dupont, E., Zhang, T., Tilke, P., Liang, L. & Bailey, W. (2018), ‘Generating realistic geology conditioned on physical measurements with generative adversarial networks’, arXiv preprint arXiv:1802.03065 .
 Emerick & Reynolds (2012) Emerick, A. A. & Reynolds, A. C. (2012), ‘History matching timelapse seismic data using the ensemble Kalman filter with multiple data assimilations’, Computational Geosciences 16(3), 639–659.
 Emerick & Reynolds (2013) Emerick, A. A. & Reynolds, A. C. (2013), ‘Ensemble smoother with multiple data assimilation’, Computers & Geosciences 55, 3–15.
 England et al. (1987) England, W. A., Mackenzie, A. S., Mann, D. M. & Quigley, T. M. (1987), ‘The movement and entrapment of petroleum fluids in the subsurface’, Journal of the Geological Society 144(2), 327–347.
 Evensen (1994) Evensen, G. (1994), ‘Sequential data assimilation with a nonlinear quasigeostrophic model using monte carlo methods to forecast error statistics’, Journal of Geophysical Research: Oceans 99(C5), 10143–10162.
 Evensen (2003) Evensen, G. (2003), ‘The ensemble Kalman filter: Theoretical formulation and practical implementation’, Ocean Dynamics 53(4), 343–367.
 Farmer (2002) Farmer, C. L. (2002), ‘Upscaling: a review’, International Journal for Numerical Methods in Fluids 40(12), 63–78.
 Fetter (2018) Fetter, C. W. (2018), Applied hydrogeology, Waveland Press.
 Garven (1985) Garven, G. (1985), ‘The role of regional fluid flow in the genesis of the Pine Point deposit, Western Canada sedimentary basin’, Economic Geology 80(2), 307–324.
 Gauthier (2014) Gauthier, J. (2014), ‘Conditional generative adversarial nets for convolutional face generation’, Class Project for Stanford CS231N: Convolutional Neural Networks for Visual Recognition .
 Gerritsen & Durlofsky (2005) Gerritsen, M. G. & Durlofsky, L. J. (2005), ‘Modeling fluid flow in oil reservoirs’, Annu. Rev. Fluid Mech. 37, 211–238.
 Goodfellow et al. (2016) Goodfellow, I., Bengio, Y. & Courville, A. (2016), Deep learning, MIT Press.
 Goodfellow et al. (2014) Goodfellow, I., PougetAbadie, J., Mirza, M., Xu, B., WardeFarley, D., Ozair, S., Courville, A. & Bengio, Y. (2014), Generative adversarial nets, in ‘Advances in Neural Information Processing Systems’, pp. 2672–2680.
 Green (1995) Green, P. J. (1995), ‘Reversible jump markov chain monte carlo computation and bayesian model determination’, Biometrika 82(4), 711–732.
 Gu & Oliver (2005) Gu, Y. & Oliver, D. S. (2005), ‘History matching of the PUNQS3 reservoir model using the ensemble Kalman filter’, SPE Journal 10(02), 217–224.
 Gulrajani et al. (2017) Gulrajani, I., Ahmed, F., Arjovsky, M., Dumoulin, V. & Courville, A. (2017), ‘Improved training of Wasserstein GANs’, arXiv preprint arXiv:1704.00028 .
 Heusel et al. (2017) Heusel, M., Ramsauer, H., Unterthiner, T., Nessler, B. & Hochreiter, S. (2017), GANs trained by a two timescale update rule converge to a local Nash equilibrium, in ‘Advances in Neural Information Processing Systems’, pp. 6626–6637.
 Higgins et al. (2017) Higgins, I., Matthey, L., Pal, A., Burgess, C., Glorot, X., Botvinick, M., Mohamed, S. & Lerchner, A. (2017), betaVAE: Learning basic visual concepts with a constrained variational framework, in ‘International Conference on Learning Representations’.
 Hinton et al. (2006) Hinton, G. E., Osindero, S. & Teh, Y.W. (2006), ‘A fast learning algorithm for deep belief nets’, Neural computation 18(7), 1527–1554.
 Hinton & Salakhutdinov (2006) Hinton, G. E. & Salakhutdinov, R. R. (2006), ‘Reducing the dimensionality of data with neural networks’, science 313(5786), 504–507.
 Holloway (2005) Holloway, S. (2005), ‘Underground sequestration of carbon dioxide—a viable greenhouse gas mitigation option’, Energy 30(1112), 2318–2333.
 Jetchev et al. (2016) Jetchev, N., Bergmann, U. & Vollgraf, R. (2016), ‘Texture synthesis with spatial generative adversarial networks’, arXiv preprint arXiv:1611.08207 .
 Kingma & Ba (2014) Kingma, D. & Ba, J. (2014), ‘Adam: A Method for Stochastic Optimization’, arXiv preprint arXiv:1412.6980 .
 Kingma & Dhariwal (2018) Kingma, D. P. & Dhariwal, P. (2018), Glow: Generative flow with invertible 1x1 convolutions, in ‘Advances in Neural Information Processing Systems’, pp. 10236–10245.
 Kingma & Welling (2013) Kingma, D. P. & Welling, M. (2013), ‘Autoencoding variational Bayes’, arXiv preprint arXiv:1312.6114 .
 Kravaris & Seinfeld (1985) Kravaris, C. & Seinfeld, J. H. (1985), ‘Identification of parameters in distributed parameter systems by regularization’, SIAM Journal on Control and Optimization 23(2), 217–241.
 Krogstad et al. (2015) Krogstad, S., Lie, K.A., Møyner, O., Nilsen, H. M., Raynaud, X., Skaflestad, B. & Others (2015), MRSTAD–an opensource framework for rapid prototyping and evaluation of reservoir simulation problems, in ‘SPE Reservoir Simulation Symposium, SPE 173317MS’, Society of Petroleum Engineers.
 Lackner (2003) Lackner, K. S. (2003), ‘A guide to CO2 sequestration’, Science 300(5626), 1677–1678.
 Laloy et al. (2018) Laloy, E., Hérault, R., Jacques, D. & Linde, N. (2018), ‘Trainingimage based geostatistical inversion using a spatial generative adversarial neural network’, Water Resources Research 54(1), 381–406.
 Laloy et al. (2017) Laloy, E., Hérault, R., Lee, J., Jacques, D. & Linde, N. (2017), ‘Inversion using a new lowdimensional representation of complex binary geological media based on a deep neural network’, Advances in Water Resources 110, 387 – 405.
 Laloy & Vrugt (2012) Laloy, E. & Vrugt, J. A. (2012), ‘Highdimensional posterior exploration of hydrologic models using multipletry DREAM (ZS) and highperformance computing’, Water Resources Research 48(1).
 LeCun et al. (2015) LeCun, Y., Bengio, Y. & Hinton, G. (2015), ‘Deep learning’, Nature 521(7553), 436.
 Leo et al. (1986) Leo, T.y., Kravaria, C., Seinfeld, J. H. et al. (1986), ‘History matching by spline approximation and regularization in singlephase areal reservoirs’, SPE Reservoir Engineering 1(05), 521–534.
 Li et al. (2018) Li, H., Schwab, J., Antholzer, S. & Haltmeier, M. (2018), ‘NETT: Solving inverse problems with deep neural networks’, arXiv preprint arXiv:1803.00092 .
 Lie (2019) Lie, K.A. (2019), An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User Guide for the MATLAB Reservoir Simulation Toolbox (MRST), Cambridge University Press.
 Lorentzen et al. (2005) Lorentzen, R. J., Naevdal, G., Valles, B., Berg, A. & Grimstad, A.A. (2005), Analysis of the ensemble Kalman filter for estimation of permeability and porosity in reservoir models, in ‘SPE Annual Technical Conference and Exhibition, SPE SPE96375MS’, Society of Petroleum Engineers.
 Lunz et al. (2018) Lunz, S., Schoenlieb, C. & Öktem, O. (2018), Adversarial regularizers in inverse problems, in ‘Advances in Neural Information Processing Systems’, pp. 8516–8525.
 Mardani et al. (2019) Mardani, M., Gong, E., Cheng, J. Y., Vasanawala, S. S., Zaharchuk, G., Xing, L. & Pauly, J. M. (2019), ‘Deep generative adversarial neural networks for compressive sensing MRI’, IEEE transactions on medical imaging 38(1), 167–179.
 Martens et al. (2012) Martens, S., Kempka, T., Liebscher, A., Lüth, S., Möller, F., Myrttinen, A., Norden, B., SchmidtHattenberger, C., Zimmer, M., Kühn, M. & Others (2012), ‘Europe’s longestoperating onshore CO2 storage site at Ketzin, Germany: a progress report after three years of injection’, Environmental Earth Sciences 67(2), 323–334.
 Miyato & Koyama (2018) Miyato, T. & Koyama, M. (2018), ‘cGANs with projection discriminator’, arXiv preprint arXiv:1802.05637 .
 Mosser et al. (2018a) Mosser, L., Dubrule, O. & Blunt, M. (2018a), ‘Conditioning of threedimensional generative adversarial networks for pore and reservoirscale models’, arXiv preprint arXiv:1802.05622 .
 Mosser et al. (2017) Mosser, L., Dubrule, O. & Blunt, M. J. (2017), ‘Reconstruction of threedimensional porous media using generative adversarial neural networks’, Phys. Rev. E 96(4), 43309.
 Mosser et al. (2018b) Mosser, L., Dubrule, O. & Blunt, M. J. (2018b), ‘Stochastic seismic waveform inversion using generative adversarial networks as a geological prior’, arXiv preprint arXiv:1806.03720 .
 Muskat (1938) Muskat, M. (1938), ‘The flow of homogeneous fluids through porous media’, Soil Science 46(2), 169.
 Muskat (1949) Muskat, M. (1949), Physical principles of oil production, McGrawHill Book Co., New York.
 Oliver & Chen (2011) Oliver, D. S. & Chen, Y. (2011), ‘Recent progress on reservoir history matching: a review’, Computational Geosciences 15(1), 185–221.
 Petzka et al. (2017) Petzka, H., Fischer, A. & Lukovnicov, D. (2017), ‘On the regularization of Wasserstein GANs’, arXiv preprint arXiv:1709.08894 .
 Plessix (2006) Plessix, R.E. (2006), ‘A review of the adjointstate method for computing the gradient of a functional with geophysical applications’, Geophysical Journal International 167(2), 495–503.
 Raftery et al. (2006) Raftery, A. E., Newton, M. A., Satagopan, J. M. & Krivitsky, P. N. (2006), ‘Estimating the integrated likelihood via posterior simulation using the harmonic mean identity’.
 Reichle et al. (2002) Reichle, R. H., McLaughlin, D. B. & Entekhabi, D. (2002), ‘Hydrologic data assimilation with the ensemble Kalman filter’, Monthly Weather Review 130(1), 103–114.
 Richardson (2018) Richardson, A. (2018), ‘Generative adversarial networks for model order reduction in seismic fullwaveform inversion’, arXiv preprint arXiv:1806.00828 .
 Roberts & Stramer (2002) Roberts, G. O. & Stramer, O. (2002), ‘Langevin diffusions and MetropolisHastings algorithms’, Methodology and computing in applied probability 4(4), 337–357.
 Roggero & Hu (1998) Roggero, F. & Hu, L. Y. (1998), Gradual deformation of continuous geostatistical models for history matching, in ‘SPE Annual Technical Conference and Exhibition, SPE49004MS’, Society of Petroleum Engineers.
 Rumelhart et al. (1988) Rumelhart, D. E., Hinton, G. E. & Williams, R. J. (1988), ‘Learning representations by backpropagating errors’, Cognitive Modeling 5(3), 1.
 Salimans et al. (2016) Salimans, T., Goodfellow, I., Zaremba, W., Cheung, V., Radford, A. & Chen, X. (2016), Improved techniques for training GANs, in ‘Advances in neural information processing systems’, pp. 2234–2242.
 Shah & Hegde (2018) Shah, V. & Hegde, C. (2018), Solving linear inverse problems using GAN priors: An algorithm with provable guarantees, in ‘2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)’, IEEE, pp. 4609–4613.
 Shao et al. (2018) Shao, H., Kumar, A. & Thomas Fletcher, P. (2018), The Riemannian geometry of deep generative models, in ‘Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops’, pp. 315–323.
 Sibson et al. (1975) Sibson, R. H., Moore, J. M. M. & Rankin, A. H. (1975), ‘Seismic pumping—a hydrothermal fluid transport mechanism’, Journal of the Geological Society 131(6), 653–659.
 Suwartadi et al. (2012) Suwartadi, E., Krogstad, S. & Foss, B. (2012), ‘Nonlinear output constraints handling for production optimization of oil reservoirs’, Computational Geosciences 16(2), 499–517.
 Tarantola (2005) Tarantola, A. (2005), Inverse problem theory and methods for model parameter estimation, Vol. 89, SIAM.
 Tikhonov & Arsenin (1977) Tikhonov, A. N. & Arsenin, V. I. A. (1977), Solutions of illposed problems, Scripta Series in Mathematics, Winston.
 Vo & Durlofsky (2014) Vo, H. X. & Durlofsky, L. J. (2014), ‘A new differentiable parameterization based on principal component analysis for the lowdimensional representation of complex geological models’, Mathematical Geosciences 46(7), 775–813.
 White (2016) White, T. (2016), ‘Sampling generative networks’, arXiv preprint arXiv:1609.04468 .
 Wu et al. (2002) Wu, X.H., Efendiev, Y. & Hou, T. Y. (2002), ‘Analysis of upscaling absolute permeability’, Discrete and Continuous Dynamical Systems Series B 2(2), 185–204.
 Yeh et al. (2016) Yeh, R. A., Chen, C., Yian Lim, T., Schwing, A. G., HasegawaJohnson, M. & Do, M. N. (2016), ‘Semantic image inpainting with deep generative models’, arXiv preprint arXiv:1607.07539 .
 Yeh & Tauxe (1971) Yeh, W. W.G. & Tauxe, G. W. (1971), ‘Quasilinearization and the identification of aquifer parameters’, Water Resources Research 7(2), 375–381.
 Yong et al. (1992) Yong, R. N., Mohamed, A.M. O., Warkentin, B. P. & Others (1992), Principles of contaminant transport in soils., Elsevier Science Publishers.
 Zhang et al. (2019) Zhang, R., Li, C., Zhang, J., Chen, C. & Wilson, A. G. (2019), ‘Cyclical stochastic gradient MCMC for Bayesian deep learning’, arXiv preprint arXiv:1902.03932 .
Appendix A Generative Model Architecture and Quality Control


Permeability transform parameters: ,
Porosity transform parameters: ,
Notation for convolutional layers:
LayerType(Number of filters), k=kernel size, s=stride, p=padding. BN=BatchNorm, PS=PixelShuffle
Parameter  Value  Unit  

Normalised Water Saturation  
Normalised Oil Saturation  
Water, Oil Viscosity  
Water, Oil Density  
Brooks Corey Exponents  
Water, Oil Compressibility  
Connate, Init., Crit. Water Saturation  
Residual, Init., Crit. Oil Saturation  
Water Injection Rate  
Bottomhole Pressure  
Reference Pressure 