Hydrocarbon reservoirs are often located in spatially complex and uncertain geologic environments, where the associated costs of drilling wells for exploration and development are notoriously high. These costs may be reduced with an optimized well-placement strategy based on real-time geologic information, known as well geosteering. To effectively place the well in an updated geomodel and support well geosteering decisions in real time, we apply an iterative inversion approach based on the Levenberg-Marquardt form of the ensemble randomized maximum likelihood method. The method estimates geomodel properties together with their uncertainties by reducing the statistical misfit between the measurements acquired with well-logging tools and the predicted measurements from numerical simulations. Analyses of synthetic cases indicate that the method’s reliability and computational speed depend on the distance from the logging tool to formation boundaries, the contrast of model properties, and the thickness of formation layers. Our method delivers reliable estimates of model properties with only 40 ensemble members and 2–10 iterations; hence, it is approximately 10–125 times faster than Metropolis-Hastings Monte Carlo, which we use as a baseline condition given its proven track record. Likewise, our method is amenable to parallelization to further reduce computational times. Implementation of the method with a synthetic example inspired by a historical well geosteering operation yields accurate formation evaluation and verifies its accurate and reliable performance under complex geologic conditions.

Recently, well geosteering has emerged as a method to optimally position wells in real time and significantly increase future production (Kullawan et al., 2014). Accurate estimation of uncertainty in geomodels while geosteering improves decisions and yields optimal well placement (Hermanrud et al., 2019). Practical well geosteering workflows were proposed by Chen et al. (2015), Kullawan et al. (2018), and Alyaev et al. (2019). All of these workflows have two main requirements: quantification of the uncertainties in model properties and a computationally efficient inversion algorithm. However, traditional deterministic and gradient-based inversion methods do not provide uncertainty quantification. In addition, correctness of the updated properties depends on the initial guess and regularization. Alternatively, methods based on Markov chain Monte Carlo (MCMC), such as Metropolis-Hastings, avoid these problems completely by iteratively reducing the error between the available measurements and the corresponding forward-simulated measurements for different geomodel realizations. Such methods deliver an optimal inverse solution and simultaneously estimate the uncertainties in the estimated modeled properties. However, such comprehensive brute-force approaches require the simulation of typically hundreds of thousands of geomodels (Metropolis et al., 1953; Bottero et al., 2016), which are computationally expensive. This condition is far from optimal for real-time analysis while geosteering.

To mitigate this problem, several research groups have proposed effective methods, such as the Hamiltonian Monte Carlo method (Shen et al., 2017) and the Bayesian inference approach (Deng et al., 2019), to reduce the sampling space without reducing the accuracy of the final solution. These approaches can be applied in a geosteering time-constrained setting, but they do not quantify uncertainties of bed-boundary locations. Ensemble-based methods, including the ensemble Kalman filter (Evensen, 1994; Aanonsen et al., 2009), are faster alternatives to MCMC, which also estimate multiple properties along with their uncertainties. Chen et al. (2015) and Luo et al. (2015) implement ensemble-based methods for joint estimation of bed-boundary locations and petrophysical properties in a geosteering setting, but their study was limited to synthetic electromagnetic (EM) tools.

In this study, we implement an iterative ensemble-based method, the Levenberg-Marquardt ensemble randomized maximum likelihood (LM-EnRML), and verify that this method is reliable for real-time inversion and uncertainty quantification of different geophysical borehole measurements (i.e., well logs). We accomplish this by estimating bed-boundary locations and petrophysical properties concomitantly with their uncertainties, across thinly layered formations exhibiting large property contrasts. In the presented examples, we use shallow-sensing nuclear density and (extra) deep-sensing borehole EM measurements. Shallow-sensing nuclear density logs are used to estimate the fine-scale properties of the beds close to the well, whereas deep-sensing EM measurements provide information about formation boundaries and properties at a coarser scale farther away from the well.

In the following section, we describe our proposed ensemble-based method. Then, using a series of synthetic examples, we verify the reliability, robustness, and computational speed of the proposed method. In one example, rock density is an uncertain model property, whereas, in the remaining examples, layer resistivity and bed-boundary locations are uncertain. Finally, we consider a formation model inspired by a historic well geosteering operation in the Goliat field in the Barents Sea, where resistivity and depths of bed boundaries are uncertain. The posterior uncertainty of the estimated model properties is quantified in all cases. We compare our results with the MCMC method (Metropolis et al., 1953), given that the latter is a well-established method with high accuracy, albeit slow.

An ensemble-based method is used to estimate the posterior mean of the uncertain model properties (i.e., formation boundaries and petrophysical properties) and to quantify their uncertainty. Here, an ensemble of geomodels, consisting of samples drawn from a Gaussian prior distribution, are conditioned to noisy measurement using Bayes’ rule, assuming that the measurement error is Gaussian. An approximate posterior distribution is obtained by solving the following minimization problem:
minxjR{12(xjxjprior)TΣx1(xjxjprior)+12(g(xj)dj)TΣd1(g(xj)dj)},
(1)
where the subscript j is the index of the ensemble member, xj is a geomodel ensemble member from the posterior, xjprior is a geomodel ensemble member from the prior, Σx is the covariance matrix for the prior ensemble, dj is a realization of the noisy measurements, and Σd is the covariance matrix for the Gaussian noise in the observed measurements. The function g(.) is the forward simulator and represents the relationship between the model properties xj and the numerically predicted measurements. When g(x) is nonlinear, equation 1 must be solved iteratively. To this end, we use the Levenberg-Marquardt algorithm.

Our predicted measurements are nuclear density measurements from the simulator developed by Mendoza et al. (2010) and (extra) deep-sensing EM measurements from a deep neural network (DNN) method (Alyaev et al., 2021) that approximates the results of an EM simulator. The DNN is trained by using data sets from a high-fidelity physics-based simulator with software provided by the tool vendor (Sviridov et al., 2014).

Algorithm

Our algorithm (Algorithm 1) is based on the method introduced by Chen and Oliver (2013). Its two main parts are the initialization and the update procedure.

Initialization

To minimize equation 1, our algorithm starts in line 2 by setting the initial guess x0,j equal to the prior ensemble xjprior. The distribution of the prior sample is Gaussian with mean xprior, covariance Σx, and Ne realizations. The initial predicted measurements d0,j are the output of the forward function g(x0,j). Afterward, in line 3, the ensemble of observed measurements dj is generated by adding realizations of Gaussian noise ϵj to the observed logging-while-drilling measurements dtrue,j. Then, in line 4, we calculate the initial data misfit and use this to set the initial value of the Levenberg-Marquardt damping factor λ (line 6).

Update

The update of the model properties from minimization of equation 1 using the Levenberg-Marquardt method at the lth iteration is given by
δxl,j=((1+λl)Σx1+Gl1TΣd1Gl1)1×[Σx1(xl1,jxjprior)+Gl1TΣd1(g(xl1,j)dj)],
(2)
where Gl1 is the sensitivity of measurements to model properties, which is calculated as the differentiation of function g(xl1) with respect to model properties xl1. By neglecting the terms containing the property mismatch and applying the Sherman-Woodbury-Morrison matrix inversion formulas (Sherman and Morrison, 1950), in equation 2 we obtain
δxl,j=ΣxGl1T[(1+λl)Σd+Gl1ΣxGl1T]1(g(xl1,j)dj).
(3)
Following Gu and Oliver (2007), we approximate Gl1 from the ensemble as
Gl1=Σd1/2ΔDl1(ΔXl1)1,
(4)
where ΔXl1 is the centered and scaled ensemble of model properties at the (l1)th iteration, and ΔDl1 is the centered and scaled ensemble of predicted measurements at the (l1)th iteration. The term ΔXl1 is calculated in line 12 as
ΔXl1=[x(l1,1),,x(l1,Ne)](INe1NeAAT)/Ne1,
(5)
where [x(l1,1),,x(l1,Ne)] is the Ne members of the ensemble of model properties, INe is the Ne × Ne identity matrix, and ARNe is the column vector of ones. Thereafter, in line 13, ΔDl1 is calculated as
ΔDl1=[g(x(l1,1)),,g(x(l1,Ne))](INe1NeAAT)/Ne1,
(6)
where [g(x(l1,1)),,g(x(l1,Ne))] is the ensemble of predicted measurement, i.e., the result of applying the forward simulator g(.) to the ensemble of model properties at the (l1)th iteration. Following Chen and Oliver (2013), we estimate Σx in the Hessian term of equation 3 from the ensemble using
Σx=ΔXl1ΔXl1T.
(7)
Note that equation 7 is a positive semidefinite matrix, and it is different in each iteration.
To stabilize the method, in line 14, we perform truncated singular-value decomposition (SVD) of ΔDl1 as
ΔDl1=UpSpVpT,
(8)
where the subscript p denotes the truncation. In this numerical study, we sort the singular values in descending order and retain the first p values such that the sum of the first p values corresponds to 99% of the sum of all singular values. The remaining singular values, with their corresponding singular vectors, are discarded. Hence, we keep 99% of the information (Anderson et al., 1999). By inserting equations 57 into equation 3 and simplifying, we obtain
δxl,j=ΔXl1ΔDl1T[(1+λl)Id+ΔDl1ΔDl1T]1Σd1/2(g(xl1,j)dj).
(9)
Insertion of the SVD from equation 8 into equation 9 gives
δxl,j=ΔXl1VpSp[(1+λl)Ip+Sp2]1UpTΣd1/2(g(xl1,j)dj),
(10)
which is the update equation in line 16. Equation 10 is repeated until the stopping criteria (c) are reached (see lines 11–21). Here, we stop the iteration when the relative improvement in data misfit is below a predefined threshold. To improve the convergence behavior, the first term of equation 2 is modified with a multiplier λ, which is calculated first in line 6 and is used in line 16. For each iteration that provides a reduction of the objective function, we set λl+1=λl/10 (line 24).

The output of the update procedure is an ensemble of models that approximately sample the posterior distribution of the model properties. The estimated model properties with their corresponding uncertainty are obtained from the mean and standard deviation of the a posterior ensemble, respectively.

In this section, we verify that the proposed ensemble-based method is reliable, robust, and computationally efficient for interpreting shallow-sensing density and (extra) deep-sensing EM well logs while quantifying uncertainty. First, we construct synthetic examples in which our method is applied to estimate layer bulk densities (using density logs), and resistivities and bed-boundary locations (using EM logs). We then consider a case inspired by the Goliat field; the method uses (extra) deep EM logs to estimate layer resistivity and bed-boundary locations. We evaluate the computational efficiency of our proposed ensemble method by performing the same example with the Metropolis-Hastings Monte Carlo method (Metropolis et al., 1953), which is the most commonly used MCMC method.

Bulk density estimation

In this example, we interpret density logs acquired across thinly laminated formations by constructing a synthetic case (see Figure 1) with layer thicknesses alternating between 0.5 and 5.5 m. Each layer has a uniform density of 2.5 (thin layers: thicknesses are less than 2 m) or 2.7  g/cm3 (thick layers: thicknesses are between 4 and 5.5 m).

The prior ensemble is a realization of density in each layer, and it is modeled with a Gaussian distribution. We run two cases, one with an ensemble of 50 realizations and another with 100 realizations. The estimated density is given by the mean of the posterior ensemble, whose standard deviation expresses the estimate’s uncertainty. Density measurements are generated from the numerical simulation of the synthetic model. We construct the measurement ensemble by sampling a Gaussian with standard deviation equal to 1.5% of the available measurements. In this example, the well path is vertical, whereas well logs are acquired at each meter along the depth of the well.

Figure 2 shows the posterior ensemble realizations (case with 50 realizations), their mean (i.e., the estimated density from LM-EnRML), the density estimated by the MCMC method, and the true value for each layer. The figure shows that, for each layer, true density lies within the posterior ensemble and is close to its mean.

Our estimates also agree well with the MCMC method, which used 10,000 sampling steps, but only three iterations of equation 10 were sufficient to achieve convergence for our LM-EnRML method. Increasing the number of realizations from 50 to 100 does not yield a significant improvement; hence, 50 realizations are considered sufficient. This yields a total of 150 forward runs compared with the 10,000 of MCMC; the proposed method is 67 times computationally less expensive than MCMC.

Shoulder-bed effects, which refer to the effect of adjacent beds on well logs, are stronger across thin beds, especially when the thickness of a layer falls below the tool’s resolution (Torres-Verdin et al., 2010; Masoudi et al., 2017). In our density simulator, shoulder-bed corrections are not applied to well logs acquired across thin beds. Therefore, the uncertainty of identifying formation properties increases across thinly bedded formations. This effect also is visible in Figure 2.

Resistivity estimation

In this example, we verify the speed and accuracy of the proposed method to estimate layer resistivities and bed-boundary locations of a geomodel with horizontally homogeneous and isotropic layers. Deep-sensing EM measurements are generated from the simulation of the synthetic geomodel. The measurement ensemble has a standard deviation equal to 1%–3% of the observed measurements. Furthermore, the well path is defined by its inclination, and well logs are acquired at each meter along the measured depth (MD) of the well.

Using three synthetic conditions, we study the sensitivity of the method to three properties: (1) distance between logging tool and layer boundaries, (2) property contrast across layer boundaries, and (3) layer thicknesses. The three synthetic cases are as follows:

  1. variable distance between logging tool and layer boundaries

  2. formation layers with high-resistivity contrast (3 versus 50 ohm-m)

  3. layers with variable thickness (0.7–10 m).

Case 1: Variable distance

In Figure 3, the tool is moving toward a high-resistivity target layer, passing through a low-resistivity layer above it. Distance from the tool to the target layer’s upper boundary decreases from 18 (at true vertical depth [TVD] 1522 m) to 1 m (TVD 1539 m). The objective is to calculate the distance from the tool to the bed boundary and quantify its uncertainty. The ensemble used for this example has 100 realizations of the depth of the bed boundary with a prior standard deviation of 2.5 m.

Figure 4 shows the posterior distribution of the estimated upper and lower boundaries of the target layer versus the distance from the logging tool. Standard deviations of each posterior distribution are shown in Figure 5; as expected, uncertainty increases when the logging tool is farther away from bed boundaries.

From Figure 5, it is clear that the uncertainty of the position of upper and lower boundaries decreases as the tool comes closer. Results indicate that the estimation of distance to the nearest (i.e., upper boundary) is accurate because the mean of the posterior realizations is always close to the true value. This behavior agrees with the study by Larsen et al. (2019), who report similar bed-boundary detection results.

Case 2: Resistivity contrast

We now consider the geomodel in Figure 6. It contains six horizontal layers with a high contrast in resistivity among them, alternating between 3 and 50 ohm-m. The first five layers from the top are drilled with a high angle (90°–80°), then the well path becomes nearly horizontal such that the bottom layer is not penetrated by the well trajectory.

We use an ensemble of 40 realizations and perform five iterations, giving a total of 200 forward runs. The plot in Figure 7 displays the resulting prior and posterior ensembles. Only in the first layer does the ensemble provide a good estimate within the 2% range of measurement noise; results are less reliable in other layers, especially those with higher resistivity. The estimated properties in layers with higher resistivity also exhibit a visibly larger ensemble spread, which is highest in the last layer located farthest away from the tool. The reason for such a high uncertainty is that, for the same EM logging tool configuration, the measured EM signal in the resistive formation cases is weaker compared with the conductive formation cases. This condition results in a reduced EM tool’s sensitivity to the high formation resistivity.

Figure 8 compares the ratio of standard deviations of prior and posterior ensembles, and it shows that the uncertainty decreases after the update. Layers with higher resistivity exhibit consistently higher standard deviation in comparison with other layers. The lowest standard deviation is found in the fifth layer and is due to the high number of measurements in this layer because it is in this layer that the well path becomes nearly horizontal (see Figure 6). The estimated resistivity in the bottom layer, farthest from the logging tool, exhibits visibly higher standard deviation.

Case 3: Variable layer thickness

To verify the method’s applicability across thin layers, we consider a synthetic case with different layer thicknesses, ranging from 0.7 to 10 m. To focus only on the effects of layer thickness, resistivity is set to 50 ohm-m in all layers. The well path inclination is 80°. We use 100 realizations, in which the mean resistivity of the prior model is 45 ohm-m and the measurement error is 3%, applied as standard deviation to the observed measurements.

Results are plotted in Figure 9, in which we observe that thin layers give rise to an estimate with a noticeably higher uncertainty. This behavior is evident from the plot of standard deviations versus layer thickness of Figure 10. Such a general trend can be explained by the low number of data points acquired within thin layers and the associated shoulder-bed effects.

Simplified field example: Goliat field paradigm

In this section, we test our method on a case inspired by a real geosteering operation in the Goliat field, located in the Barents Sea, as described by Larsen et al. (2015). The geomodel, shown in Figure 11, covers a small section without faults before the well lands in the drilling target layer in the Upper Kobbe Formation. The formation includes thin layers and high-resistivity contrasts. Compared to the actual formation, the earth model and the well path have been rotated 5° to make the layering horizontal. The well path is from MDs 1000 to 1105 m, sampling a total of 106 positions along the well path. Well inclination increases with depth from 80° to 85°.

Unknown model properties are layer resistivities and boundary depths, and they are estimated simultaneously. We construct the measurement ensemble by sampling a Gaussian with standard deviation equal to 1% of the observation value. The prior ensemble has 100 realizations, with standard deviations of 3% for resistivities and 0.25% for bed boundaries (which, for the depths considered in Figure 11, is approximately 2.5–3.2 m).

Resistivity estimation

Figure 12 shows the standard deviation of the resistivities of the prior and posterior ensemble. In general, uncertainties decrease after the update, although not significantly for the two bottom layers, which are located farther away from the tool. This behavior is expected because the sensitivity of the EM tool decreases with distance to the layers.

The posterior resistivity distribution for each layer is shown in Figure 13, normalized to the true value. True resistivity values of layers 4 and 6 lie inside the interquartile range of the posterior ensemble, with layer 3 barely outside.

Layer 5, which is the first layer below the tool, has an estimated resistivity between the prior estimate and the true value, and the lowest standard deviation of all layers. This result is consistent with layer 5 having the lowest true resistivity. As observed in case 2, low resistivities give rise to lower uncertainties because the EM response from conductive formations is relatively stronger. Conversely, the second layer’s high standard deviation is due to its high resistivity, in addition to it being a thin layer.

Figure 14 compares the estimated resistivity from the ensemble method with the resistivity obtained using the Metropolis-Hastings Monte Carlo method (Metropolis et al., 1953) and with the true resistivity. Estimated resistivities using the Metropolis-Hastings Monte Carlo exhibit higher accuracy for the first two top layers. For layers 3 and 4, both methods show comparable accuracy because there are enough data from the measurement sampling points. Neither of the methods are reliable in estimating the resistivity of the bottom layer due to lack of data.

The ensemble-based method requires 3–5 iterations, yielding 300–500 runs of the forward simulators. By comparison, Metropolis-Hastings Monte Carlo requires at least 10,000 runs of the forward simulator to minimize the data misfit and estimate the resistivity with comparable accuracy.

Bed-boundary estimation

The ensemble method can simultaneously estimate multiple unknowns. In this example, we estimate bed-boundary locations along with layer resistivities. The relative standard deviation of bed boundaries in the prior ensemble is 0.25%, equivalent to approximately 2.5–3.2 m.

The posterior standard deviation, shown in Figure 15, significantly decreases compared to the prior one for boundaries 2–5. Boundaries 1 and 6–8, which are located farther away from the tool are, as expected, more uncertain.

In this paper, we proposed an ensemble-based method for petrophysical inversion based on the iterative LM-EnRML and verified it on several petrophysical inverse problems. It can simultaneously estimate bed-boundary locations and layer petrophysical properties from noisy measurements, in addition to providing uncertainty estimates at no extra computational cost. Results indicate that the method reduces uncertainties in the depths of layer boundaries and petrophysical properties within the depth of investigation of well-logging tools. Estimation results are less precise in geologic models that include thin layers and rock formations with highly resistive layers, and in layers farther away from the logging tool.

The main advantages of the proposed method are (1) the number of forward simulations required is at least 20 times lower than for a standard Metropolis-Hastings MCMC method while yielding similar results and (2) the method is amenable to parallelization, making it especially attractive for computationally expensive applications when performance time is critical. Such applications include real-time ensemble-based interpretation while geosteering, training of machine-learning models, real-time testing of geologic hypotheses against the measurements, and other applications in which a large number of forward simulations is needed within a limited time.

Finally, we emphasize that our implementation enables the possibility of joint inversion of multiple well logs in the same framework, enabling consistent, unbiased, and effective interpretation of multiple types of while-drilling measurements for real-time geosteering decision support.

The authors thank S. Tveit, senior researcher at NORCE, for valuable discussion, and the reviewer for the constructive feedback. The authors also gratefully acknowledge the support from the research project “Geosteering for IOR” (NFR-Petromaks2 project no. 268122), which is funded by the Research Council of Norway, Aker BP, Equinor, Vår Energi, and Baker Hughes Norway. Part of this work was carried out within the financial support from the Center for Research-based Innovation DigiWells: Digital Well Center for Value Creation, Competitiveness and Minimum Environmental Footprint (NFR SFI project no. 309589, DigiWells.no). The center is a cooperation of NORCE Norwegian Research Centre, the University of Stavanger, the Norwegian University of Science and Technology (NTNU), and the University of Bergen, and funded by the Research Council of Norway, Aker BP, ConocoPhillips, Equinor, Lundin, TotalEnergies, and Wintershall Dea.

Data associated with this research are available and can be obtained by contacting the corresponding author.

graphic
Nazanin Jahani received a Ph.D. (2015) in mechanical engineering from the Norwegian University of Science and Technology (NTNU) in Trondheim. From 2010 to 2017, she contributed to the development of a near wellbore reservoir simulator for Petrell AS in Trondheim. Since 2017, she has been a research scientist at NORCE. She is a project member for the project “Geosteering for IOR,” working on statistical methods and applying them for real-time petrophysical interpretation.

graphic
Joaquín Ambía received a Ph.D. (2011) in computational biophysics from the University of Houston. From 2011 to 2014, he worked on the Rosetta software for computational modeling of structural biology systems. Since 2014, he has been a research associate and software developer in the Formation Evaluation Joint Industry Research Consortium of the University of Texas at Austin. He contributes to the 3D UTAPWeLS well-logging simulation and earth-model integration software. He is responsible for incorporating newly developed algorithms and developing new algorithms and methods.

graphic
Sergey Alyaev received an M.S. and a Ph.D. in applied mathematics from the University of Bergen. He is a senior researcher in applied mathematics at NORCE with interests including forward and inverse modeling, data-driven methods, and decision theory. The main focuses of his research are the real-time applications of these methods and integrated workflows in the field of drilling and geosteering. From 2018, he took responsibility for leading the project “Geosteering for IOR” where the team is developing advanced quantitative and statistical methods and applying them for real-time decision making.

graphic
Kristian Fossum received an M.S. in petroleum technology and a Ph.D. in applied mathematics, both from the University of Bergen. He is a researcher in the data assimilation and optimization group at NORCE. His research interests include ensemble-based data assimilation, inverse problems, computational statistics, data-driven methods, optimization, and automated decision support.

graphic
Erich Suter received a B.S. in information technology from Gjøvik University (now NTNU), an M.S. in applied mathematics from the University of Oslo, and a Ph.D. in computational geoscience from the University of Stavanger. He is a senior research scientist at NORCE. He is skilled within applied mathematics and computer science, petroleum geology, geometric modeling, and software design and development. His research interests are within real-time geologic interpretation and earth model management, for applications such as decision support under uncertainty while geosteering and automatic construction/updating of geomodels.

graphic
Carlos Torres-Verdín received a Ph.D. (1991) in engineering geoscience from the University of California at Berkeley. During 1991–1997, he held the position of research scientist with Schlumberger-Doll Research. From 1997 to 1999, he was reservoir specialist and technology champion with YPF (Buenos Aires, Argentina). Since 1999, he has been affiliated with the Hildebrand Department of Petroleum and Geosystems Engineering of the University of Texas at Austin, where he is currently full professor. He is the founder and director of the Research Consortium on Formation Evaluation at the University of Texas at Austin. He is the recipient of the Cockrell School of Engineering’s 2016–2017 Lockheed Martin Aeronautics Company Award for excellence in engineering teaching, the 2014 Gold Medal for Technical Achievement from the SPWLA, the 2008 Formation Evaluation Award from the SPE, and the 2006 Distinguished Technical Achievement Award from the SPWLA. He is a distinguished member of the SPE, honorary member of the SEG, and received the Conrad Schlumberger Award from the EAGE.

Freely available online through the SEG open-access option.