DeepFlow History Matching in the Space of Deep Generative Models

History Matching in the Space of Deep Generative Models

Lukas Mosser
Department of Earth Science and Engineering
Imperial College London
&Olivier Dubrule
Department of Earth Science and Engineering
Imperial College London
\ANDMartin J. Blunt
Department of Earth Science and Engineering
Imperial College London

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 ill-posed 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 pre-trained on geostatistical object-based 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 two-phase incompressible Darcy formulation. We invert for the underlying reservoir properties by first modeling property distributions using the pre-trained 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 rock-type 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 non-linear forward problem.


History Matching in the Space of Deep Generative Models

  Lukas Mosser Department of Earth Science and Engineering Imperial College London Olivier Dubrule Department of Earth Science and Engineering Imperial College London Martin J. Blunt Department of Earth Science and Engineering Imperial College London


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), carbon-capture 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 non-aqueous 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 multiphase-fluid 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 long-term 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 ill-posed inverse problem.

For an inverse problem to be well-posed it must satisfy three criteria (Tikhonov & Arsenin, 1977)

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

  2. Uniqueness: the obtained solution has to be unique.

  3. 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 ill-posed. It has been shown that inverse problems with spatially distributed parameters common in groundwater and reservoir flow problems are often ill-posed 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 well-defined. 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 B-splines 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 optimisation-based 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 history-matching 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 non-linear problems.

Recent developments in deep learning (LeCun et al., 2015; Goodfellow et al., 2016) have motivated the use of neural-network-based 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 (ES-MDA) (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 Markov-Chain Monte-Carlo methods (Laloy & Vrugt, 2012). Chan & Elsheikh (2018a) used generative adversarial networks trained using patch-based kernel discrepancy measures to synthesise geological models from training images. Realisations conditional to spatial observations of the underlying model parameters can be obtained by gradient-based 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).

Adjoint-state 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 Metropolis-adjusted 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 (ill-posed) 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 non-linear 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 two-phase immiscible Darcy flow equations without gravity and under isothermal conditions using a finite-difference approach and obtain gradients of the mismatch between observed and simulated data with respect to the gridblock permeability and porosity using the adjoint-state method.

Starting from random initial locations within the underlying Gaussian latent space controlling the GAN’s output, we use gradient-descent 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 two-dimensional synthetic case of a channelised reservoir system we show that gradient-based inversion using GANs as a model parameterisation can be achieved, honoring the observed physical quantities of the non-linear forward problem as well as borehole measurements of the underlying petrophysical quantities. The trained models and code are available as open-source software111

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 bottom-hole pressure. The forward problem is therefore defined as solving the two-phase (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 .


We consider all fluid parameters of the forward problem, the oil and water viscosity , density , compressibility , as well as the relative-permeability behaviour, represented by a Brooks-Corey 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 single-phase permeability , and porosity .


We assume a binary set of reservoir rock-types where corresponds to a low permeability shale and indicates the high permeability sandstone of a river-channel.

Given the non-linear forward operator this allows us to define the forward problem as


We represent the distribution of possible model parameters by a (deep) generative latent-variable model


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:


and from Bayes theorem


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


In the case of GANs only individual samples can be obtained and therefore estimation of the posterior is restricted to point-estimates. We derive a regularisation scheme for the history matching problem from finding point estimates of the maxima of the posterior distribution (MAP)


where the squared-error 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 cross-entropy between the observed facies indicator and the facies probability of each grid-block 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 squared-error functional that measures the difference between the observed pressures and rates at the well locations averaged over the total duration of the observed data


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


and the standard deviation of the pressure measurements


used to generate the set of observed pressure and rate data.

Using the adjoint-state 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


consisting of the flow loss (Eq. 9), the well-data 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 non-linear forward problem as well as the likelihood function (Eq. 9) have been previously implemented in the open-source MATLAB/GNU-Octave reservoir simulation framework MRST (Lie, 2019). An overview of the end-to-end coupling between the latent variables, the deep generative model , and the numerical forward model is shown in figure 1.

Figure 1: Overview of a single forward and backward-pass through the combined generative network and forward problem of two-phase Darcy flow . A set of model parameters is generated from a forward-pass through the generative model . The generative model creates continuous spatial properties such as the permeability shown here (top). The forward problem is solved where water is injected in the injection well (blue) and oil and water are produced from the production well (red). The error functional between observed parameters and modelled observations is evaluated. A spatial distributed gradient (bottom) with respect to the underlying properties is computed using the adjoint state method. The obtained gradients are then backpropagated by traditional neural-network backpropagation to compute a gradient with respect to the latent variables. The error functional is minimised by performing gradient descent using the obtained gradient.

To evaluate the proposed approach we create synthetic two-dimensional vertical cross-sections of a stacked river-channel system using an object-based approach. Channel-bodies are represented by half-circles with centers following uniform distributions in space. The full object-based modeling workflow is detailed in Mosser et al. (2018b). For the case of reservoir history matching (Eq. 3) each object-based model is associated with the three model parameters following equation 2 (Fig. 2).

Figure 2: Object-based model used to generate reference observed data. Each object-based model has three associated properties: the rock-type (shale 0 - sandstone 1) indicator (a), the grid block permeability (b), and porosity (c). The model was obtained from a test-set of object-based model realisations that the generator network was not trained on.

We train a deep generative adversarial network (GAN) (Goodfellow et al., 2014) on a set of object-based realisations to learn the generative model (Eq. 4) prior to obtaining MAP samples (Eq. 8) by optimising equation 12. Specifically, we train a so-called Wasserstein-GAN (Arjovsky et al., 2017; Gulrajani et al., 2017), that minimises the Wasserstein distance between generated and real probability distributions. To stabilise training of the Wasserstein-GAN we use the Lipschitz penalty introduced by Petzka et al. (2017). A detailed overview of the training procedure for the deep-generative 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 100-dimensional noise-sample from a standardised multi-Gaussian distribution to a spatial distribution of the model parameters (Eq. 2) with 128 pixels in the x- and 64 pixels in the z-direction


where the three output features correspond to the grid-block rock-type 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 rock-type can be interpreted as a rock-type probability. The grid-block permeability and porosity values are obtained by transforming the generated rock-type 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 finite-difference approximation of the forward model is computed on the regular grid given by the pixel-based output of the generative model and hence each pixel-based property corresponds to the grid-block property used in the numerical evaluation of the forward problem.

We have selected a reference realisation from a test-set of object-based 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. 10-11). 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 gradient-based 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

Figure 3: (a) Overview of the dynamic flow behaviour, as well as average and standard deviation of the gridblock rock-types for samples obtained from the unconditional prior distribution i.e. samples drawn at random from the generative model. (b) Model realisations after optimising only the gridblock rock-type at wells, (c) optimising only the dynamic flow behaviour, and (d) optimising both well rock-type and flow behaviour. Red contours on the mean and standard deviation of rock-type distributions reflect the reference model of object-based channel-bodies shown in figure 2.

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. 3-a).

  • 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. 3-b).

  • Scenario 3: We performed inversion by considering only the flow-based loss due to the mismatch between simulated and observed dynamic two-phase flow data, as well as the prior loss on the latent variables (Fig. 3-c).

  • Scenario 4: All three losses were combined: the flow loss , the well-data loss , and the prior loss , corresponding to the total loss , as outlined in equation 12 (Fig. 3-d).

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 4: Histograms of the total loss functional (a, c) and the grid-block rock-type accuracy at the wells (b, d). Models () optimised on the observed flow data only (a-b, red) and optimised on both flow and well rock type (c-d, red) are compared to models () sampled unconditionally from the prior distribution (a-d, black). In both cases, we observe that optimising the flow-based objective leads to realisations that match the observed data closer than sampling unconditionally from the prior. The accuracy of the well rock type for models optimised only on the dynamic flow data follow the prior distribution of models, whereas the models obtained by optimising flow and well-data loss functions have high accuracy at the well locations.

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 rock-type probability maps. The threshold for rock-type 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. 3-a, right) shows a large spread in the distribution of the observed oil-water 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 rock-type indicator function , are shown as cross-sections of the mean and standard deviation of the models (Scenario 1, Fig. 3-a, 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 grid-block rock-type indicator function at the wells only (Scenario 2) where the binary cross-entropy 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. 3-b, right). Inspecting the mean and standard deviation maps of the rock-type 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 grid-block well data can be observed comparing the outline of the mean rock-type 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. 3-c, right). The mean and standard deviation maps for the obtained distribution of the rock-type 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.

Figure 5: Comparison of the number of models that have a connected cluster of channel-bodies that connects the injection and production wells in each individual sample obtained in the inversion process for the four scenarios outlined in section 2. After Boolean dilation, nearly all models that were optimised with respect to the observed dynamic flow data form a connected cluster between the two wells, whereas half of the samples obtained from honoring only the well data remained disconnected.

Combining the flow loss , the well data loss and prior loss as the total loss (Eq. 12) (Scenario 4, Fig. 3-d) 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 grid-block data at the production and injection well . This is indicated by the low standard deviation map of the inverted models (Scenario 4, Fig. 3-d, left).

Figure 6: Optimisation trajectory for a selected inverted sample based on flow, well-data and prior losses. The flow loss shows a highly non-linear optimisation problem with many local minima leading up to the smallest error being achieved after 214 iterations. The well data accuracy at the production and injection wells quickly reaches values greater than 90% accuracy. Cross-marks indicate the loss values of the intermediate realisations presented in figure 8.

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 well-data loss is optimised in addition to the flow-loss (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 forward-model 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 total-loss 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).

Figure 7: Comparison of the observed dynamic pressure and rate data (red) obtained for a single inversion case. The observed dynamic data from the initial starting model (green) is optimised by honoring flow and well data for 500 iterations. The intermediate pressure and rate data of the 500 iteration steps are colored according to their total loss .

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. 3-c, 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. 3-b, 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. 13-14). We evaluate for all four scenarios whether any of the clusters of river-channel 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 finite-difference approximation (two-point flux approximation) used to solve the forward-problem 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 high-permeable channel-bodies 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 latent-variables (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 rock-type indicator function at the well locations (red), and the prior loss evolution. A highly non-linear 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.

Figure 8: Evolution of the spatial distribution of the rock-type indicator for the selected model.

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 a-b). The solution found by the optimisation strategy (Eq. 8) with the lowest error value (Fig. 7 a-c, black) closely matches the observed dynamic production data and the rock-type indicator function at the wells.

The intermediate spatial distributions of the rock-type indicator for each of the obtained local minima in the optimisation process are shown in figure 8. Initially the channel-bodies were completely disconnected (Fig. 8 a). As the optimisation process progresses more channel-bodies were introduced gradually (Fig. 8 b-e) and form a spanning cluster that matches the grid-block data at the wells (Fig. 8 f).

4 Discussion

We have shown that solutions of the ill-posed 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 grid-block 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 multi-modal, where each of the obtained realisations lie at a minimum of the total loss i.e. near a local maximum of the posterior.

(a) Total Loss and well facies accuracy for 100 realisations interpolated between three obtained maxima of the posterior .
(b) Flow Loss , well loss , and prior loss of the interpolated latent variables between three maxima obtained by optimising the total loss (Scenario 4).

While the obtained samples match the flow and well data it is important to note that in a high-dimensional 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 history-matched 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 river-channel 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 three-dimensional multi-rock-type models i.e. more than two rock-types 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 net-to-gross. 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 class-conditional variables to GANs, so-called cGANs, result in improved training and distributional representation behaviour compared to their unconditional counter-parts (Gauthier, 2014). This may speak for training a single general reservoir model generator compared to creating a library of standalone pre-trained 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, mode-collapse, and representation of spatial statistics as highlighted in Mosser et al. (2017) and Laloy et al. (2018), where two-point 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) pre-trained networks and the features observed in typical reservoir architectures. Kernel-based 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 gradient-based approach to solving the ill-posed inverse problem is feasible, other strategies such as using Markov-Chain Monte-Carlo 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 Kalman-Filters (Gu & Oliver, 2005; Aanonsen et al., 2009) or Ensemble Smoothers (ES-MDA) (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 batch-size of one, it is possible to perform inversion using a batch-size 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 latent-variables and hence for a set of model parameters . By computing the numerical solution forward problem and the associated adjoint for each set of model-parameters 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 ill-posed inverse problems opens up interesting avenues with regards to the theoretical basis of the ill-posed 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 latent-space of deep generative models and Shao et al. (2018) have shown that the manifolds obtained from VAEs are non-linear with near-zero 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 ill-posed inverse problems. Nevertheless, GANs due to their challenging training, evaluation and complex latent-to-parameter-space relationship may not be the optimal reparameterisation choice in the case of ill-posed 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 flow-based 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 ill-posed inverse problems.

5 Conclusions

We have presented an application of deep generative models in the context of the ill-posed history matching inverse problem. Based on a two-dimensional synthetic cross-section of a river-channel system we have trained a GAN to represent the prior distribution of the subsurface properties , permeability , porosity and the rock-type indicator . We find solutions to the two-phase slightly-compressible Darcy flow problem of oil and water systems commonly used to describe the transient flow behaviour in hydrocarbon reservoirs. By incorporating a finite-difference-based numerical simulator with adjoint-state capabilities (Krogstad et al., 2015) in an end-to-end differentiable framework, we perform inversion using gradient-based optimisation of the mismatch between observed dynamic data (Eq. 9) and grid-block-scale well data (Fig. 3), while honoring the prior distribution of the latent variables (Eq. 12). By using a momentum-accelerated first-order gradient descent scheme (Kingma & Ba, 2014), the method converges to a solution despite the highly non-linear and non-convex 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 latent-variable generative models and their theoretical benefits for finding solutions to ill-posed inverse problems.


O. Dubrule would like to thank Total S.A. for seconding him as visiting professor at Imperial College London.


  • (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 ill-posed inverse problems using iterative deep neural networks’, Inverse Problems 33(12), 124007.
  • Adler & Öktem (2018) Adler, J. & Öktem, O. (2018), ‘Learned primal-dual 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., Maier-Hein, 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 pore-scale 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), ‘Pore-scale 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 Learning-Volume 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 beta-VAE’, 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 XVI-16th European Conference on the Mathematics of Oil Recovery’.
  • Chadwick & Noy (2010) Chadwick, R. A. & Noy, D. J. (2010), History-matching flow simulations and time-lapse 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), ‘Exemplar-based 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 steady-state groundwater flow’, Journal of Hydrology 31(3-4), 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 141510-MS’, 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. & Fei-Fei, L. (2009), ‘Imagenet: A large-scale hierarchical image database’.
  • Dinh et al. (2016) Dinh, L., Sohl-Dickstein, 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 time-lapse 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 quasi-geostrophic 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(1-2), 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., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, 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 PUNQ-S3 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 time-scale 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), beta-VAE: 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(11-12), 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), ‘Auto-encoding 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), MRST-AD–an open-source framework for rapid prototyping and evaluation of reservoir simulation problems, in ‘SPE Reservoir Simulation Symposium, SPE 173317-MS’, 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), ‘Training-image 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 low-dimensional 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), ‘High-dimensional posterior exploration of hydrologic models using multiple-try DREAM (ZS) and high-performance 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 single-phase 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 SPE-96375-MS’, 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., Schmidt-Hattenberger, C., Zimmer, M., Kühn, M. & Others (2012), ‘Europe’s longest-operating on-shore 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 three-dimensional generative adversarial networks for pore and reservoir-scale models’, arXiv preprint arXiv:1802.05622 .
  • Mosser et al. (2017) Mosser, L., Dubrule, O. & Blunt, M. J. (2017), ‘Reconstruction of three-dimensional 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, McGraw-Hill 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 adjoint-state 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 full-waveform inversion’, arXiv preprint arXiv:1806.00828 .
  • Roberts & Stramer (2002) Roberts, G. O. & Stramer, O. (2002), ‘Langevin diffusions and Metropolis-Hastings 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, SPE-49004-MS’, Society of Petroleum Engineers.
  • Rumelhart et al. (1988) Rumelhart, D. E., Hinton, G. E. & Williams, R. J. (1988), ‘Learning representations by back-propagating 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 ill-posed 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 low-dimensional 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., Hasegawa-Johnson, 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

Latent Variables
Conv2D(512)k3s1p1, BN, ReLU, PSx2
Conv2D(256)k3s1p1, BN, ReLU, PSx2
Conv2D(128)k3s1p1, BN, ReLU, PSx2
Conv2D(64)k3s1p1, BN, ReLU, PSx2
Conv2D(64)k3s1p1, BN, ReLU, PSx2
Conv2D(64)k3s1p1, BN, ReLU, PSx2
= permeability
= porosity
(a) Multi-Property Generator
Model Parameters
Conv2D(64)k5s2p2, ReLU
Conv2D(64)k5s2p1, ReLU
Conv2D(128)k3s2p1, ReLU
Conv2D(256)k3s2p1, ReLU
Conv2D(512)k3s2p1, ReLU
Conv2D(512)k3s2p1, ReLU
Conv2D(1)k3s1p1, ReLU
(b) Discriminator
Table 1: Generator and discriminator network architectures used to represent the distribution of model parameters used in the inversion process (Eq. 2). Binary indicators of geological facies and derived permeability are represented by a bi-variate Gaussian distribution with a hyperbolic tangent activation function used to create the binary output distribution. We use a linear transformation layer to renormalise the output of the GAN from the interval of the tanh activation to a known range of permeability and porosity values.
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
Table 2: Fluid and simulation parameters for oil and water phases used in the numerical solution of the two-phase flow forward problem.
Figure 9: Distribution of the number of optimisation iterations required to reach a total loss (Eq. 9) of less than or equal to in the case of matching flow and well rock type data (Scenario 4). This threshold corresponds to the lower end of the total-loss distribution of the prior distribution (Fig. 4, c). Optimisation converges on average within the first 30 iterations.
Figure 10: Comparison of the pdfs of samples of the test set and samples obtained from the generative network trained to represent the river-channel body system. The reservoir rock-type (a) and permeability (b) show a near binary distribution whilst porosity follows a bi-modal pdf. The pdfs of the generated samples match the test set data closely.

Appendix B Inverted Samples

Figure 11: Unconditional samples () obtained from sampling the prior distribution of the generative model (Sec. 3, Scenario 1)
Figure 12: Samples obtained by minimising the sum of the well and prior losses (Sec. 3, Scenario 2)
Figure 13: Samples obtained by minimising the sum of the flow and prior losses (Sec. 3, Scenario 3)
Figure 14: One hundred samples obtained by minimising the sum of the flow, well, and prior losses (Sec. 3, Scenario 4)
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description