Verification of a real-time ensemble-based method for updating earth model based on GAN

The complexity of geomodelling workflows is a limiting factor for quantifying and updating uncertainty in real-time during drilling. We propose Generative Adversarial Networks (GANs) for parametrization and generation of geomodels, combined with Ensemble Randomized Maximum Likelihood (EnRML) for rapid updating of subsurface uncertainty. This real-time ensemble method combined with a highly non-linear model arising from neural-network modeling sequences might produce inaccurate and/or biased posterior solutions. This paper illustrates the predictive ability of EnRML on several examples where we assimilate local extra-deep electromagnetic logs. Statistical verification with MCMC confirms that the proposed workflow can produce reliable results required for geosteering wells.


Introduction
The process of drilling wells for hydrocarbon production represents a major cost in petroleum reservoir development. However, drilling of new wells is necessary to increase the total oil recovery. To maximize the value for each drilled well it is necessary to optimize the placement of the well within the reservoir structure. An optimally placed well will mobilize more of the petroleum resources, and reduce the need for injected water -reducing the environmental impact of oil production.
To place a well in its optimal position, operators apply geosteering. Here, the well trajectory is adjusted while drilling in response to real-time measurement of the geology surrounding the drill bit. The value of geosteering has been well documented in the literature [1,2,3].
The main objective with geosteering is to utilize the information in the measurements to make optimal decisions. Hence, geosteering can be seen as a sequential decision process under uncertainty and should be treated in a probabilistic framework [4]. Recently, a workflow based on the Ensemble Kalman Filter (EnKF) [5] has been employed to condition the geological model on measurements acquired while drilling [6,7]. In the EnKF, the uncertainty is represented by an ensemble of equiprobable realizations. This workflow has then been combined with a global optimization method and applied as a Decision Support System (DSS) [8].
The DSS framework provides high quality decisions on synthetic cases, and outperforms most of geosteering experts in a controlled experiment [9]. However, practical challenges should be addressed for it to be applicable to real operations [8]. This includes modeling of modern commercial tool to process real measurements as well as real-time earth model that can handle realistic geological complexity. The forward deep neural network (FDNN) trained on synthetic data for extra-deep electromagnetic measurements [10] enabled the real-time ensemble-based update of layered models in 1.5D [11], and 3D [12]. Moreover, [13] showed that the model errors present in the FDNN approximation can be alleviated during the ensemble based inversion for the layered case.
Fossum et. al [14] proposed a new modeling sequence which combines the FDNN with a generative adversarial network (GAN) to produce complex geological realizations in real-time to aid geosteering, see Figure 1. The premise of the GAN is that it allows to represent the earth model by a Gaussian distribution, where all produced realizations also maintain geological realism. This allows using EnKF-like methods to update not only continuous properties but also complex geological structures, which is required for geosteering. However, the workflow implementation in [14] converged only with a little starting uncertainty. [15] improved the results and demonstrated visually-convincing structural ahead-of-bit prediction on a selected example. It is known, however, that the approximate ensemble based methods, such as EnKF and its derivatives, can be sensitive to non-linearities present in complex modeling sequences and thus predictions may be biased.
In this paper we present a robust and improved implementation of the framework presented in [14] that is able to account for an appropriate starting uncertainty. Further, we aim to test the convergence properties of the iterative ensemble randomised maximum likelihood (EnRML) method when updating the GAN-based geomodels with FDNN approximation of measurements -denoted the GAN-FDNN modeling sequence. The EnRML probabilistic output is compared to a gold standard Markov Chain Monte-Carlo (MCMC) solution for the same problem using various metrics. The numerical examples demonstrate that the EnRML -applied to the GAN-FDNN modeling sequence -generates posterior samples with excellent predictive capabilities that are good approximations to the true posterior solution.
To construct a reference earth model we generate realizations of a fluvial geological environment using a commercial software. These realizations are then sub-sampled to form a training dataset for the offline training of a Generative Adversarial Network (GAN). The GAN is then used, online, to generate plausible geological realizations from a low-dimensional Gaussian input vector. The complete earth modeling loop is described in Section 2. For modeling the extradeep EM measurements we use a forward deep neural network (FDNN) trained on a dataset generated using a commercial simulator (Section 3). In Section 4 we discuss the exact and the approximate data assimilation (DA) methods. The two numerical experiments, designed to test the applicability of our proposed method, are derived, and the numerical results are presented in Section 5. Finally, we summarize and conclude the paper in Section 6.

Earth modeling using GAN
GANs are a class of unsupervised machine learning methods which can learn to generate new formatted data with the same statistics as the training set. Motivated by successful applications of GANs for modeling channelized structures for reservoir-simulation workflows [16,17,18,19,20], we use a GAN for efficient earth modeling.
The GAN consists of two deep neural networks (DNNs): a generator and a discriminator. The generator takes a random Gaussian low dimensional vector as input and generates a realization of formatted data: geological realization. The discriminator takes the formatted data and gives a probability of it being 'real', i.e., belonging to the training set. During training the DNNs contest each other in a min-max game. They are trained simultaneously. On each training step the generator creates (fake) geological realizations from the random vectors. Fake geological realizations are combined with random samples of the real earth model and are fed to the discriminator. The loss function for the generator is proportional to number of 'fakes' correctly identified by the discriminator. The loss function for the discriminator is proportional to the total misjudged data samples. In our study we use an adapted Wasserstein GAN [21] with hierarchical deep convolutional networks [22] for the generator and the discriminator, see [21] for implementation details.
For geosteering we want to reproduce likely geological realizations of facies and porosity distributions on a 2D vertical geological section along the well to identify the oil-bearing sands ahead of bit. For training of the GAN we use a large (compared to the area of prediction) reference earth model, which should provide a realistic test case for the present study in terms of scale and actual geological features and properties. The reference earth model is constructed using a commercial software that models a synthetic structural framework, a facies model setup derived from outcrop analogue data, and synthetic petrophysical properties of individual facies derived from published literature. The resulting model measures 4000m x 1000m x 200m (xyz) with cell dimensions set to 10x10x0.5 m, yielding a regular grid of size 400x100x400, see Figure 2.
The constructed facies model represents a low net/gross fluvial depositional system. It was chosen since it provides complex 3D architectures comprising a limited number of facies, which form contrasting geometries, see Figure 2. Input numbers for statistical generation of facies and geometries are derived from a well-documented outcrop of the Cretaceous lower Williams Fork Formation (Mesa Verde Group) at Coal Canyon, Colorado, USA [23,24,25].
Key parameters of the facies model setup are listed in Table 1. The model is not intended as a rendering of the outcrop itself and is consequently simplified compared to descriptions of the original outcrop [26,27,28,24]. The model contains three facies: Background/shale, Channels and Crevasse splays. The probability distribution of channel width in the model is adapted to include "narrow channel bodies", and stacking of channels accounts for multistory channels which comprise more than 80% of the observed channel bodies. The flow direction of the channel system is set towards 45 ± 10 degrees. No trends were used to condition the spatial distribution of channels. The details of this synthetic model are also described in [15]. The geological realization is parameterized by a vector of 60 independent parameters. For each 60-dimensional vector, the generator outputs a 64x64 grid with three values in each grid block. For a grid block (with dimensions 10.0m along-well and 0.5m thickness) the three values, 'channels', represent the probability of the grid-block belonging to the respective facies class: Background/Channel/Crevasse. Our generator is also predicting porosity/resistivity distribution within the geo-bodies, but in this study only the facies classes are used.
For training, the original 3D earth model is sampled as 64x64 2D images with three channels. The facies index from the training set is converted into one-hot three-dimensional vector. That is, the vector represents the probability of facies: the value of the true index is set to one and other channels to zero. During evaluation, the resistivity of the facies with the highest probability is applied.

Forward DNN model of extra-deep EM logs
To maintain real-time performance of a data assimilation workflow the forward model should be fast and support batch, preferably parallel execution. Proprietary forward models provided by measurement instrument vendors provide the most accurate results, but they are often not sufficiently fast, and not always optimized for batch execution. In [10], the authors developed a DNN approximation of such a forward model [29], which we abbreviate FDNN.
The model approximates the output of the ultra-deep electromagnetic wellbore logging instrument. The instrument is configured to transmit four shallow and nine pairs of deep directional measurements, and has sensitivity to boundaries up to 30 meters to the side from the well bore. We emphasize that the tool provides information around, but not ahead of the drilling position. An illustration of the deep measurements depth of detection is provided in Figure 11.
The input to the FDNN model is a layered geological media with up to three boundaries above and below the measurement instrument as well as the resistivity values of all seven layers. In this study we assume that the layer resistivity is isotropic and that the well is aligned with the horizontal axis.
We produce one synthetic set of measurements for every horizontal position (one per column of cells) of the gridded model which we 'drill' through. We choose the most probable facies for each computational cell within the considered column and use the corresponding resistivity value (same as in [15]): We find the boundaries between layers composed of pixels with equal resistivities and use the boundaries and the layers' resistivities as the input to the forward model.

Data assimilation during geosteering
The DSS for geosteering [8] uses the data assimilation loop (see Figure 1, left) to condition the earth model to measurements made while drilling. The fundamental idea is that if a poorly known earth model can be made consistent with measurements in the statistical sense, it will contain non-biased forecasts and, hence, provide a better basis for decisions (see Figure 1, right).
In this paper, the emphasis is placed on the data assimilation part of the DSS. Specifically, we investigate real-time data assimilation with the EnRML method utilizing a modeling sequence based on the two neural networks described above: a GAN-generator for complex earth modeling and a FDNN models for the synthetic extra-deep EM logs. The EnRML method is an ensemble-based iterative ensemble smoother which has received a lot of attention for history matching subsurface multi-phase flow problem, see [30] and references therein. The method uses an ensemble approximation to the sensitivity matrix, and provides a fast and approximate solution to the Bayesian problem. The method can only be shown to converge for Gaussian posterior distributions. However, the method is known to also sample accurately from moderately non-Gaussian posterior distributions. To assess the statistical convergence of the EnRML in this context we compare the method to samples generated by a Markov Chain Monte Carlo (MCMC) algorithm. Since the MCMC, when properly converged, generates independent and identically distributed samples from the Bayesian distribution, several metrics can be used to assess the statistical error of the EnRML method. The EnRML and the MCMC method is described in detail in the rest of the section.

EnRML
The EnRML [31] has recently become one of the most successful methods for automatic history matching of petroleum reservoirs. The EnRML is derived by minimization of an objective function using an ensemble approximation of the sensitivity matrix. Since a wide range of methods can be applied for the minimization, the EnRML can be formulated in many different ways. In this study we utilize the approximate form of the Levenberg-Marquardt method, introduced in [32].
Based on the Bayes' theory, the objective function to be minimized is Here, d * obs is the noisy observed data, g ( m) is the modelling sequence depending on the parameter m, and m prior is a sample from the prior distribution of the parameters (this is the rough plus smooth sampling approach given in [33, chap.10]). Iteration number i of the Levenberg-Marquardt method is given as where λ i is the Levenberg-Marquardt multiplier, G is the sensitivity of data to the parameters, and ∼ N 0, C d is a realization of the measurement observation noise.
In the ensemble framework, we approximate C m and G using the ensemble. To this end we defineG where N denotes the ensemble size, and C sc is a diagonal matrix for scaling the data, typically containing the measurement variance on the diagonal. We get the approximate version of the Levenberg-Marquardt update equation by inserting ensemble approximations of G and C m , neglecting the updates from the model mismatch term, substituting the prior precision matrix C −1 m withC −1 mi , and rewriting the equation using the Sherman-Woodbury-Morrison matrix inversion formula [34] gives the following update equation The update equation is simplified, and made computationally more stable, by inserting the truncated singular value decomposition of ∆d where the subscript p indicates the number of retained singular values, when they are ordered after descending value. In this work, we define p such that the cumulative sum of the p first singular values equals 99% of the cumulative sum of all the singular values. Further, we substitute C D with the ensemble approximationC DC where Inserted into (8) gives where Z and ζ are the eigenvectors and eigenvalues of After each application of (8) we assess convergence and we continue iterations until the scheme is is converged. Here, we consider the method to be converged when the relative difference in the data misfit (the first term in (1)) is below a given threshold or when the maximum number of iteration is reached. In the numerical examples, the threshold is 2 × 10 −2 and the maximum number of iterations are 10.

MCMC
A reliable method for sampling from a complex posterior distribution is the MCMC technique. MCMC relates to the general framework of methods introduced in [35] and [36] for Monte Carlo (MC) integration. Firstly, one designs a Markov chain that produce samples from the desired posterior distribution. Secondly, one utilize these samples for MC integration. In this section, the adaptive Metropolis-Hastings method -the method utilized in the numerical study -is introduced. For more information on MCMC we refer the reader to [37], and references therein.
Suppose we want samples from the un-normalized posterior distribution F , which is the general case with the Bayesian method where the normalizing factor often is very difficult to calculate. Assume that the current element of the chain is m, and one proposes a move to m * . The proposal is sampled from the proposal distribution q ( m * | m). The move is performed with probability b ( m, m * ) = min (1, r ( m, m * )) (13) where the Hastings ratio is defined as This is the basis for the Metropolis-Hastings method, and it can be shown that the method generates samples from the posterior distribution F . The Metropolis-Hastings algorithm requires a choice of proposal distribution, and some distributions work better than others. The optimal would be to draw proposals directly from the posterior F . However, this is not possible since we cannot sample from this distribution. Since the MCMC converges for any proposal distribution that fulfills some general conditions, one idea is to gradually adapt the proposal distribution using previous samples from the chain. This adaptive approach ensures a gradually better proposal distribution as the chain evolves. To this end we select the following mixture distribution as our proposal distribution HereC m is the empirical covariance matrix calculated utilizing all the preceding iterations of the Markov Chain, Q m is some fixed non-singular matrix and 0 < β < 1. Note that β = 1 untilC m is well defined. Efficient on-line updating of C m is achieved by the recursion given in [38]. This sampling method was applied in [39,40]. It is well known that the MCMC requires a certain burn-in period since the initial samples are not from the posterior distribution. Hence, it is necessary to monitor the convergence of the method. In this work, convergence is monitored by assessing the maximum root statistic of the multivariate potential scale reduction factor [41].

Numerical Experiments
The numerical experiments investigate how the GAN-FDNN modeling sequence can be applied in the data assimilation part of DSS when data assimilation using EnRML is applied for real-time uncertainty reduction. We design two synthetic experiments that focus on the reduction of uncertainty ahead of measurements and the ability of the algorithm to predict the sand channels in the unexplored part of the geomodel.
To quantify the quality of the EnRML approximation, when applied to the GAN-FDNN modeling sequence, we compare the posterior ensemble of EnRML with true samples from the posterior -acquired by the MCMC algorithm. The comparison is done by evaluating several metrics, including visual comparison of standard deviation, and mean, in every point of the domain, point-wise Kolmogorov-Smirnov two-sample test, and visual inspection of kernel density estimates of the marginal distribution of GAN-input vector m i .
We perform two numerical tests. In both tests we utilize the generative neural network, introduced in Section 2, to represent the earth model with uncertainty. Hence, our goal is to condition the poorly known 60-dimensional input vector, m, to measurements. The prior realizations of the earth model are generated by applying the generative network to parameters sampled from a multivariate Gaussian distribution, m ∼ N ˆ 0 m, C m . The distribution is slightly shifted to simulate conditioning on pre-drill information. Theˆ 0 m represents the shifted mean and is defined by the equation equation: where m 0 is the synthetic truth from [15], see Figure 6. We use uncorrelated covariance matrix with marginal variance of C mi = 1 for all parameters, similar to the GAN training. Figure 3 shows six generated earth model realizations from the prior model. From the figure, we observe that this setup provides significant variation in the earth model, which is reflected in the relatively flat mean and the standard deviation derived from the full ensemble of 500 realizations, see Figure 5. At the same time, all the realizations are consistent with the chosen channelized geological setting. We emphasize, that we use the same prior for both numerical experiments. We conduct two numerical experiments, that differ with respect to the synthetic truth. In the first experiment, the truth model from [15] is applied. In the second experiment, the synthetic truth depicted in Figure 11 is applied. Hence, for experiment 2, the prior is biased towards a wrong model (indicating erroneous pre-drill information) making the data assimilation problem harder.
The numerical study is performed in the same manner for both experiments. Firstly, we sample the true posterior with the MCMC. Here, 8 Markov chains, starting from different initial points, were run for 10 6 steps. At that point, based on assessing the multivariate potential scale reduction factor, the MCMC was found to be converged. Samples from the posterior were then extracted by removing the burn-in phase, and by thinning. For each of the 8 chains, the first half of the chain was removed, and every 100th iteration from the second half of the chain was retained, leaving 4 × 10 4 samples from the posterior distribution. Secondly, we estimate the posterior distribution using the EnRML method introduced in Section 4.1. Due to the fast simulation time, we utilized an ensemble size of N = 500, and in addition, we applied the correlation-based localization technique introduced in [42]. Finally, we assess the result from the EnRML by comparison with the samples from the MCMC. 1.

5.
6.    The resistivity of an earth model generated by GAN used as the synthetic truth for Example 1 (adapted from the numerical example in [15]). The yellow arrows show the region with measurements and their extent illustrates the maximum sensitivity range, termed depth of detection. The filled red line is the drilled well, and the dashed lines indicate the potential for geosteering.

Example 1 -Verification of convergence on an example from literature
The first numerical example tests the sampling capabilities for the EnRML on an example from the literature. The synthetic true log is generated from the true model 6. Hence for this case, the prior mean is slightly shifted towards the true model. Figure 7 shows the mean and standard deviation of the posterior resistivity model. Compared to the prior mean and standard deviation (shown in figure 5) it is clear that the uncertainty around the well is significantly reduced. Moreover, the same reduction is observed for both the EnRML and the MCMC. Apart from slightly sharper boundaries for the MCMC, the EnRML approximation to the posterior mean and standard deviation is almost indistinguishable from the true posterior mean and standard deviation.
A similar conclusion can be drawn from the plots showing the estimated point-wise probability of the sand facies, plotted in Figure 8. Around the well,  the EnRML and the MCMC are almost indistinguishable. Ahead of the bit, there are small differences. However, the predictive capability, as described in [15], is also present in the MCMC solution. Hence, this is a true feature of the posterior solution with the GAN-FDNN modeling sequence.
To evaluate the statistical distance between the samples from the EnRML and the samples from MCMC we perform a Kolmogorov-Smirnov two-sample test. This is a non-parametric test of equality for one-dimensional probability distributions. The earth model is a 2D image, and not one-dimensional. Hence, we perform the test on the marginal distribution for each cell. The P-values from the test of the H0 hypothesis of equal distributions are shown in Figure 9. The significance level of 0.05 is given by the black contour line. Hence, for all p-values higher than 0.05 one cannot differentiate between the marginal distributions and the H0 hypothesis hold.
As a final evaluation of the results, we plot a selection of marginal and bivariate elements of the input vector m. To highlight the effect of the pre-drill information, we selected two elements that was shifted (m 1 and m 52 ) and two that were not shifted (m 23 and m 37 ). In Figure 10 the kernel density estimate of the selected elements is plotted along the diagonal, the scatter plot of the pairwise elements is given in the top corner, while the contours of the 2D Kernel density estimate of the pairwise elements are given in the lower corner. The true model is given as a black line in the 1D plots and a black star for the 2D plots.
The numerical experiment shows that the EnRML can successfully approx-

Example 2 -Prediction of a sand-channel sequence
The second numerical example tests the ability of the workflow to predict the targets ahead of measurements in the case where the well is already landed into a sand channel, and when the pre-drill information, embedded in the prior model, is biased toward a wrong solution. The synthetic truth for this example with the depth of detection is shown in Figure 11. Figure 12 illustrates the posterior mean and standard deviation of the resistivity model. It is clear that conditioning to measurements resolves the sand channel ahead of the well position. Moreover, the EnRML does a reasonably good job in approximating the true posterior, despite the prior being slightly misspecified.
Similarly, the estimated point-wise probability of the sand facies, plotted in Figure 13, shows that the EnRML provides excellent predictive capabilities as the sand facies is correctly forecasted to the right of the geomodel, more than 500 meters ahead of the bit. There are slightly larger differences between the a. b.   EnRML and MCMC in this example. However, the approximate posterior is still very close to the MCMC posterior, especially around the well. The measure of statistical distance between the samples from the EnRML and the samples from MCMC indicates similar performance. The P-values from the test of the H0 hypothesis of equal distributions are shown in Figure 14. The significance level of 0.05 is given by the black contour line. Hence, for all p-values higher than 0.05 one cannot differentiate between the marginal distributions and the H0 hypothesis hold. Compared to example 1, there are more areas where the H0 hypothesis fails. This demonstrates that, for the more challenging experiment, there is a larger statistical distance between the approximate posterior from the EnRML and the true posterior.

EnRML MCMC
As a final evaluation of the results, we plot a selection of marginal and bivariate elements of the input vector m. Similar to experiment 1, we highlight the effect of the biased pre-drill information by selecting two elements that were shifted (m 14 and m 53 ) and two that were not shifted (m 20 and m 342 ). In Figure 15 the kernel density estimate of the selected elements is plotted along the diagonal, the scatter plot of the pairwise elements is given in the top corner, while the contours of the 2D Kernel density estimate of the pairwise elements are given in the lower corner. The true model is given as a black line in the 1D plots and a black star for the 2D plots. The effect of the biased prior can be observed in the MCMC results, where the marginal posterior is bi-modal.
The numerical experiment shows that the EnRML can successfully approx-imate the true posterior solution for the GAN-FDNN modeling sequence, even when the prior model is slightly biased. The numerical results show convincing similarity between the exact samples from the posterior, acquired by the MCMC, and the approximate samples, acquired by EnRML. There is however a larger discrepancy than was observed in example 1. From inspection of selected elements from m we can observe that the biased prior results in a more non-Gaussian posterior distribution. Despite this, we claim that the EnRML provides good approximations of both the posterior earth model and the posterior input vector.

Conclusions
In this paper, we have demonstrated that two essential parts, the earth model and the simulated extra-deep EM logs, of an ensemble-based DSS system can be substituted with neural networks. For the earth model, we utilize the GAN trained with images from a realistic geological setting, for the simulated logs we use a forward deep neural network (FDNN) trained using a large set of simulations from a commercial tool. The setup redistributes the computational cost from online to offline calculations, enabling complex earth models and deepsensing EM logs to be part of real-time ensemble updates.
The numerical results illustrate that the GAN-FDNN modeling sequence provides excellent probabilistic predictions ahead of drilling capturing both continuous and discrete features when conditioning to only measurements with sideways sensitivity. Moreover, the numerical results show that the computationally efficient EnRML algorithm can sample the true Bayesian posterior confirmed by the MCMC algorithm. This conclusion is valid even when the prior model is slightly biased towards a wrong solution.
The proposed approach has many beneficial factors. Firstly, a GAN provides large flexibility for defining the geological setting. Here, we consider three different facies, but one can easily imagine including features like faults and pinch-outs as well as smoothly-varying properties. Secondly, we only need to condition a few parameters with Gaussian distribution to the measurements, which is very beneficial for the ensemble-based DA approach. Thirdly, since we are utilizing a neural network model to generate the simulated log, the computational cost of simulating a single ensemble member is milliseconds. Hence, the proposed approach can utilize a large ensemble for the DA part.
The numerical experiments illustrated that the posterior has a predictive capability for both MCMC and the faster EnRML method. The future work is to integrate the DA developed in this paper with the decision framework developed in [8], allowing DSS under a much more complex geological setting. Furthermore, the method can be extended to account for model errors present in machine learning approximations in real-time [13].