Bayesian joint inversion of seismic and electromagnetic data for reservoir lithofluid facies, including geophysical and petrophysical rock properties

INTRODUCTION

The goal of mitigating risk in hydrocarbon exploration and reservoir recovery has driven the development of tools to integrate seismic data, well logs, core samples, and geologic models in reservoir characterization (Grana et al., 2012). There are various approaches for quantitative estimation of reservoir properties from seismic inversion; some invert for elastic and petrophysical properties as continuous parameters (Bosch et al., 2010), and others use facies estimation (Buland and Omre, 2003; Grana and Rossa, 2010; de Figueiredo et al., 2019a; Pendrel and Schouten, 2020; Grana et al., 2022). However, mapping fluid saturation from seismic signals faces critical limitations because the rock matrix dominates the acoustic and elastic impedances (Mavko et al., 2009). Such limitations make the interpretation highly dependent on prior knowledge (Buland et al., 2008), signal processing, and acquisition (Zerilli et al., 2018).

In contrast, the electrolytic conductivity of fluids in the connected pore spaces dominates the resistivity of sedimentary rocks (Mavko et al., 2009). Indeed, resistivity well logs are the primary tool for indirectly measuring fluid saturation throughout a borehole. In particular, the observed strong contrasts in resistivity of a siliciclastic reservoir at the oil-brine interface have led to the increasing use of electromagnetic geophysical methods for mapping fluids in marine reservoirs (Buonora et al., 2014; Alvarez et al., 2018; Miotti et al., 2018; Menezes et al., 2021).

Despite its relatively recent history, the marine controlled-source electromagnetic (CSEM) method has been successfully applied to offshore hydrocarbon exploration (Constable, 2010), often at the intermediate stage of field development: after seismic interpretation because CSEM acquisition is focused on leads and the interpretation of the acquired signals depends on prior information, and before reservoir characterization (Buonora et al., 2014), due to the limited signal resolution compared to the typical geologic scale in which reservoir engineers are interested.

Some approaches invert the stand-alone CSEM data for a resistivity model that focuses on the seismically derived zone of interest and then use that resistivity and elastic rock properties from seismic inversion to estimate petrophysical properties such as porosity, saturation, shale volume, and lithofluid facies by calibrated rock-physics modeling (Alvarez et al., 2018; Miotti et al., 2018).

However, constraining resistivity models with seismically derived regularization and prior models may need to be revised or more accurate in recovering models at the inner reservoir scale. For instance, assigning transverse resistance for the interval from CSEM-only inversions does not contribute to the seismic mapping of the oil-water contact (OWC) and requires a circular interpretation process (Andréis et al., 2018). Therefore, the seismic resolution to lithology should be associated with CSEM sensitivity to fluid saturation by simultaneous joint inversion linked in the model domain (Zhdhanov et al., 2012).

Techniques of seismic and CSEM joint inversion have been proposed for linking elastic to electrical models by spatial constraints (Hu et al., 2009), petrophysical coupling (Hoversten et al., 2006; Chen et al., 2007; Gao et al., 2012; Tu and Zhdanov, 2022), and facies categorization (Chen and Hoversten, 2012; Pendrel and Schouten, 2020). Inverting for facies links electric to elastic models by spatial discretization associated with facies-related prior means and covariances of geophysical rock properties. Petrophysical coupling can be explicitly introduced when using petrophysical parameters as prior distributions and associating rock physics with geophysical modeling in the forward function (de Figueiredo et al., 2019b) or, as in this work, implicitly when using petrophysical parameters transformed by rock-physics modeling to compose the prior distributions.

Over the past few decades, several works have been concerned with facies inversion from geophysical data, often presenting the Bayesian approach as a natural framework to solve inverse problems that combine categorical and continuous variables. The posterior probability density function (PDF) defines the target solution for the model parameters, providing the maximum posterior probability (MAP) model and related uncertainty. Given the number of cells in 3D reservoir models (Zerilli et al., 2018), the full stochastic evaluation of the posterior PDF (Chen and Hoversten, 2012; de Figueiredo et al., 2019b) might be very costly when modeling CSEM and seismic signals by numerical solvers of differential or integral equations (Fagin, 1991; Avdeev, 2005; Valente et al., 2021). Although this work uses 1D models and forward functions, we are interested in an inversion routine prepared to efficiently incorporate 3D models and functions.

We propose a Bayesian joint inversion of seismic and CSEM data to estimate models of lithofluid facies, in contrast with Chen and Hoversten (2012), by using a linear Gaussian approach in which the marginalization over geophysical rock properties is found analytically (Grana et al., 2017; de Figueiredo et al., 2019a). Although the Zoeppritz reflectivity (Mavko et al., 2009) and CSEM 1D modeling (Crepaldi et al., 2011) are nonlinear functions, they are used as forward modeling and assumed to be linear approximations around the facies-related averages to take advantage of such analytical solutions. Numerical tests validate this adaptation of the linear Gaussian approach for the given conditions.

The petrophysical prior distributions converted to geophysical rock properties by calibrated rock-physics modeling may introduce facies-related means and covariances different than directly acquired geophysical rock properties. Hence, we assemble prior distributions of geophysical and transformed petrophysical rock properties in the marginalization integrals. We also linearize the rock-physics modeling around the averages. An introduced parameter weights these distributions, providing stand-alone geophysical (de Figueiredo et al., 2019a) or petrophysical priors at the edges and joint distributions (Bosch, 2004; Bosch et al., 2010; Grana and Rossa, 2010; Chen and Hoversten, 2012; de Figueiredo et al., 2017) at the middle of its range.

The proposed linearizations using analytical Jacobians of forward functions efficiently evaluate the data-conditioned likelihood for a given facies model (posterior likelihood). However, there are $F^N$ possible configurations for F facies arranged in N layers, and computing them may be impracticable. Thus, a Markov chain Monte Carlo (MCMC) algorithm (Metropolis et al., 1953) samples models through that likelihood to provide a full posterior distribution in a reasonable run time.

As in de Figueiredo et al. (2019a) and different than in Chen and Hoversten (2012), the model layer thickness is equivalent to the seismic data sampling, which means a model discretization on the limit of signal resolution. Unlike these works, we do not impose spatial correlation between layers in the covariance matrices to focus on evaluating the inversion accuracy for each given data set. However, we set correlations between layers in the model perturbation.

A siliciclastic marine reservoir in the Campos Basin in Brazil (Zerilli et al., 2018) was selected to validate the methodology in a realistic scenario. We used two wells, one as prior information and another to base the synthetic model for simulating the geophysical data. After presenting and discussing the well’s data, rock-physics calibration, and matching between priors and the synthetic model, we apply the proposed lithofluid facies inversion to the simulated data for many prior settings and compare stand-alone with joint inversions.

To challenge the CSEM inversions, we simulate two pessimistic scenarios, one with lower resistivity contrasts for the fluid saturation and another with a resistive layer in the overburden, which is ignored during forward modeling. We also test many data sets in these simulations.

The linear Gaussian approach allows us to compute the conditional distributions of geophysical and petrophysical rock properties by sampling linear transformations over their respective means and covariances on each accepted facies model (Buland and Omre, 2003; Grana and Rossa, 2010; de Figueiredo et al., 2017). Instead, because we use nonlinear forward functions, nonlinear inversions (Tarantola, 2005; Pujol, 2007) are performed to minimize the residuals in these local inversions. Finally, our method is applied to joint geophysical data and joint prior distributions, providing profiles of posterior distributions for the rock properties in three inversion experiments, each using one of the three wells as a prior, including the testing well itself. We approach this feature as a supplementary subroutine without testing many parameterizations and detailing the results.

METHODOLOGY

This section discusses a Bayesian approach under Gaussian assumptions and linearized modeling for estimating the posterior distribution of the facies model within an interval of interest (IOI), including the proposed prior distribution that integrates geophysical and petrophysical rock properties, along with the numerical implementation of the theory using an MCMC algorithm (Metropolis et al., 1953) for facies inversion. The calculation and implementation of local inversions for rock properties also are presented as an optional subroutine.

Statistical model

The desired posterior PDF $p(\pi |d)$ of the facies model $\pi$ conditioned on the geophysical data d can be expressed by the following Bayesian inference (Tarantola, 2005):

where $p(d|\pi )$ is the PDF for any possible geophysical data conditioned on a proposed facies model, and $p(\pi )$ is the PDF for any possible facies models and may be conditioned on prior information about facies spatial distribution within the IOI.

The term $p(d)$ normalizes the joint distribution by the following summation over the entire space of possible facies configurations:

such that for the model space $\mathrm{\Omega }=\{{\pi }_1,{\pi }_2,\dots ,{\pi }_{n_f}\}$, the $n_f$ discriminated facies arranged in n layers result in $n_f^n$ possible configurations (Grana et al., 2017).

The PDF of the geophysical data conditioned on the facies model $p(d|\pi )$ is given by the product of the data likelihood $p(d|m)$ with the prior distribution $p(m|\pi )$ marginalized in the m properties (Grana et al., 2017; de Figueiredo et al., 2019a):

Assuming Gaussian data noise and introducing the notation ($N_{\mathrm{data}}$ [function, covariance]) for normal distributions (ND) (Buland and Omre, 2003), the data likelihood of equation 3 can be written as

where $C_d$ is the data noise variance matrix and $g(m)$ is the geophysical forward modeling. This work defines the geophysical properties m as compression $V_{\mathrm{p}}$ and shear $V_{\mathrm{s}}$ seismic velocities, density Dens, and resistivity Res. See Appendix A for more details about the seismic and CSEM forward modeling.

The prior PDF, $p(m|\pi )$, in equation 3 can be computed from the set of m properties directly acquired from well log or core sample measurements for each facies. Thus, under the Gaussian assumption, it can be written as

where ${\mu }_m$ is the vector of prior means and $C_m$ is the covariance matrix. Although not indicated by subscripts for simplicity, all means and covariances of rock properties are facies-related variables.

Alternatively, the prior PDF in equation 3 also can be computed from the petrophysical hard data r (porosity $\varphi$, water saturation $S_w$, and volume of shale $V_{\mathit{sh}}$) when it is transformed by calibrated rock-physics modeling $m_r=t(r)$ (see Appendix B). In this case, the prior PDF of the transformed geophysical properties conditioned on facies takes the form of the following marginal distribution:

Then, assuming PDFs in equation 6 as NDs produces

where ${\mu }_r$ is the vector of means, $C_r$ is the covariance matrix of the r prior distributions, and $C_{m_r}$ is a diagonal matrix with the squared errors of the rock-physics calibration, also assigned to each facies.

The following exponential balancing is now proposed to generalize equation 3 by including the priors of equations 5 and 7:

where the parameter $0\le \alpha \le 1$ is introduced to weight each prior distribution type under the requirements of retrieving the stand-alone priors at the edges of this range and ensuring that $\alpha$ never decreases the covariances of either distribution, avoiding overfitting prior averages. The normalization of this combined prior distribution is described in Appendix C.

Thus, substituting equation 8 in equation 3 leads to the following generalized marginalization:

Appendix C assumes that rock-physics $t_{\mathrm{Lin}}$ and geophysical $g_{\mathrm{Lin}}$ modeling are linear functions around the facies-related means ${\mu }_r$ and ${\mu }_m$ to solve the marginalization integrals 7 and 9 analytically. However, in the following, we substitute the nonlinear functions $t(r)$ and $g(m)$ in equations C-7 and C-13 because it may improve the accuracy.

Hence, inserting ${\mu }_{m_r}=t({\mu }_r)$ into equation C-7 and $g({\mu }_{\eta })$ into equation C-13 adapts the nonlinear forward functions to $p(d|\pi )$:

where C is the integrated covariance matrix:

The derivation of equations 10 and 11 can be validated by assuming the particular case $\alpha =0,1$ because the same results are obtained by directly substituting $p_m(m|\pi )$ or $p_r(m|\pi )$ in equation C-11 and then applying equation C-1.

In equation 11, the precision (inverse of the covariance) matrices related to each type of prior distribution are summed, ensuring that in the middle of the range ($\alpha =0.5$), the more accurate the information between the priors, the more significant its contribution to $p(d|\pi )$.

After setting a realistic geologic parameterization, the linearization and nonlinear adaptation are implicitly validated for the rock-physics modeling through the almost $\alpha$-independent results in the “Facies from priors” subsection and explicitly for the geophysical modeling in the “Geophysical linearization” subsection.

Implementation

In the following paragraphs, the theory described in the previous subsection is applied to estimate the posterior facies model distribution, conditioned on the observed geophysical data and accounting for the proposed prior modeling, along with local inversions for geophysical and petrophysical rock properties.

Facies from geophysical data

The analytical solution in equation 10 eliminates the need for numerical integration over the properties m (equation 3) for each facies configuration. However, the posterior distribution $p(\pi |d)$ (equation 1) can only be entirely determined by the evaluation of the product $p(d|\pi )p(\pi )$ over the entire facies model space ${\mathrm{\Omega }}^n$ (equation 2), which may be numerically unfeasible. Therefore, an MCMC algorithm draws multiple facies realizations from this high-dimensional posterior distribution (de Figueiredo et al., 2019a).

The probability for each proposed facies model $p(\pi )$ in equation 1 is assumed to be a first-order Markov chain:

where the facies probability at a given layer l depends on the adjacent facies at position $l-1$ according to the prior probabilities of facies occurrence $p({\pi }_1)$ and transitions $p({\pi }_l|{\pi }_{l-1})$. These probabilities are geologic parameters based on prior knowledge about the spatial distribution of the facies in the IOI (de Figueiredo et al., 2019b; Talarico et al., 2020).

At each iteration, a model window is randomly selected, keeping the same facies as the previous model in the top and bottom layers of the window. Then, the layers in between are perturbed using a truncated Gaussian simulation, according to the matrix P, the elements of which are the facies transition probability between adjacent layers $p({\pi }_l|{\pi }_{l-1})$ in equation 12. Because those probabilities are uniformly applied over the layer pairs, P is a square matrix with dimensions equal to the number of facies (de Figueiredo et al., 2019b). This windowing is adopted to increase the acceptance rate compared to the perturbation of the entire model.

The geophysical rock properties are computed on a logarithmic scale and then concatenated on a vector by blocks related to each layer. Seismic and CSEM data also are concatenated into a data vector. The Jacobians and covariance matrices are built by layer-related blocks consistent with the rock properties vector.

After evaluating $p(d|{\pi }^{\mathrm{new}})$ (equation 10) for the current model proposal, the acceptance rate is given by

where $p(d|{\pi }^{\mathrm{last}})$ relates to the last accepted model. Because $p(\pi )$ is already introduced in model perturbation, it is not present in the acceptance rate.

Then, according to the Metropolis algorithm (Metropolis et al., 1953; de Figueiredo et al., 2019a), a random number y with a uniform distribution in the range $0\le y\le 1$ is generated and compared with r. If $y\le r$, the current configuration ${\pi }^{\mathrm{new}}$ is accepted. Otherwise, it is rejected, retaining the last previous model (Figure 1). If accepted, the model sample is stored along with its averages ${\mu }_{\eta }$, responses $g({\mu }_{\eta })$, covariance $C_{\eta }$, and Jacobians G and T, and ${\pi }^{\mathrm{new}}$ becomes ${\pi }^{\mathrm{last}}$ on the next iteration (Figure 2).

In contrast with Chen and Hoversten (2012) and de Figueiredo et al. (2019a), no spatial correlation between layers is imposed on the covariance matrices to focus on the accuracy of each geophysical data set under different prior settings and $\alpha$ effects (equation 8).

Beyond the burn-in period, stored samples provide a posterior distribution for which the MAP is the most accepted model, and the uncertainties are extracted for each layer independently.

Local inversions

In the linear Gaussian approach, geophysical and petrophysical rock properties can be estimated as conditional distributions by sampling local inversions over the accepted facies models (Grana et al., 2017; de Figueiredo et al., 2019a):

where the means and covariances conditioned on the data and facies model for the proposed prior distribution (subscript $\eta$) are given by the following linear relations: and

However, under the assumed paradigm of using nonlinear forward functions, local inversions for rock properties achieve better adjustments through the following iterative process (Marquardt, 1963; Tarantola, 2005; Pujol, 2007; Miotti et al., 2018):

where k is the iteration counter, ${\mu }_{m|d,\pi }^0={\mu }_{\eta }$, $\lambda$ is the Marquardt parameter, and I is the identity matrix (Marquardt, 1963; Pujol, 2007). The inversion stops when model changes reach a minimum value $|{\mu }_{m|d,\pi }^{k+1}-{\mu }_{m|d,\pi }^k|<\epsilon$. The covariance of the posterior properties $C_{m|d,\pi }$ is updated with the final Jacobian $G_k$ in equation 16.

The same process is used for the petrophysical properties ${\mu }_{r|m,\pi }$, by substituting d into the final ${\mu }_{m|d,\pi }$, the Jacobian G to T, $C_{\eta }$ to $C_r$, and $C_d$ to $C_{m_r}$ and similarly to its posterior covariance $C_{r|m,\pi }$. The properties “r” are limited within the critical limits by a logarithm barrier (Frisch, 1955).

In this work, the two local inversions are done in cascade for quality control. However, the petrophysical properties could be directly inverted from the geophysical data by the composite function $g(t[r])$, changing ${\mu }_{\eta }$ and $C_{\eta }$ to ${\mu }_r$ and $C_r$, respectively, and G to GT and $C_d$ to $C_d+G{C}_{m_r}{G}^T$, respectively, in equation 17, which has the advantage of coupling the elastic properties to resistivity and the disadvantage of not computing the geophysical rock properties.

To compute the posterior distribution $p(\gamma |d)$ of the rock properties $\gamma \in \{m,r\}$ (with m in logarithm and r in linear scale according to the inversion parameterization), an M-element (i-index) linearly spaced vector is created for each property $\gamma$ and layer (with ranges that cover more than one standard deviation). The NDs $N_{{\gamma }_i}({\mu }_{\gamma j},C_{{\gamma }_j})$ are calculated with the means ${\mu }_{{\gamma }_j}$ and covariances $C_{{\gamma }_j}$ given by local inversions for each accepted facies configuration represented by the j-index, which includes repetitions in configuration acceptance. Thus, the average of those NDs over the j-index provides an M dimension posterior probability vector for each property and layer.

APPLICATION

The inversion routine is applied to a realistic geologic scenario to assess the proposed hybrid prior modeling and the robustness of the inversion in dealing with many different prior settings. Three wells of a marine reservoir were selected, one for the synthetic model and the others to provide the prior distributions.

This section presents the well data, the calibration of rock-physics parameters, and the ability of the priors to predict the reference facies from the synthetic model of geophysical rock properties. Then, facies inversions are applied to synthetic data, scanning the $\alpha$ range for each prior well and geophysical data set. Two possible challenges in CSEM applications are simulated: inversions for the lowest possible resistivity contrast and using a misguided overburden model.

The local inversions for rock properties are applied as a potential resource without scanning the prior settings or interpreting the results at the same level of detail as in the facies inversion.

Well data

Three Maastrichtian deep-water reservoir wells (confidential geographic reference) are selected for this work: WP1, WP2, and WT, where (W) refers to the well, (P) to prior, and (T) to test. Each well contains the three lithofluid facies used in this work: shale, brine sand ($S_w\ge 0.5$), and oil sand ($S_w\le 0.5$), the last two separated by the indicated OWC, all inside the previously interpreted top and bottom of the IOI (Figure 3).

Notably, the resistivity of oil-saturated sandstone presents a high anisotropy ratio (Santos and Régis, 2015). However, the conventional resistivity well logs of the vertical wells used are only sensitive to the horizontal resistivities ${\mathit{Res}}_h$. Meanwhile, the sensitivity of the isotropic CSEM modeling used (see Appendix A) is dominated by the vertical resistivity ${\mathit{Res}}_v$ (Newman et al., 2010).

Because there is no triaxial induction well log on these selected vertical wells (Abubakar et al., 2006), the ratio of electric anisotropy for each facies (${\theta }_{\pi }$) is coarsely estimated from horizontal wells and a fluid flow simulator on the reservoir (Zerilli et al., 2018): ${\theta }_1=2$, ${\theta }_2=1.5$, and ${\theta }_3=10$. The vertical resistivity profiles are then computed by applying $\theta$ as a multiplicative factor on the resistivity data from the three wells (Figures 2 and 4). The fluid flow simulator also provides the effective pressure ($P_{\mathit{eff}}$) used in rock-physics modeling (see Appendix B).

Figure 4 shows WT original and upscaled lithofluid facies (shale in gray, brine sand in yellow, and oil sand in green) with 13 layers of 5.4 m thickness each, determined by the adopted seismic time sampling (see Appendix B), along with geophysical and petrophysical properties from WT well logs, including the rescaled ${\mathit{Res}}_v$, the synthetic model, and the prior averages and standard deviations related to the upscaled facies model. The upscaling algorithm builds the synthetic rock properties model by averaging values from the original profiles inside the respective upscaled layer. A quality control uses the inference method described in the “Facies from priors” subsection to ensure the synthetic model of geophysical rock properties has the highest probability of belonging to the upscaled facies model.

Figure 5 shows crossplots for the petrophysical properties $V_{\mathit{sh}}$ versus $\varphi$, elastic $I_s$ versus acoustic $I_p$ impedances, and vertical resistivity ${\mathit{Res}}_v$ versus water saturation $S_w$, for WPs and WT. Each contour corresponds to one standard deviation of property pair distributions in multivariate form. The thinner contour lines are geophysical rock properties transformed from petrophysical properties by calibrated rock-physics modeling (see Appendix B).

Shale presents sharp and poorly correlated distributions due to the shortage of well-log samples inside the IOI (Figure 3). To better characterize shale, it is necessary to complement the priors with well data samples outside the IOI limits, possibly discriminate shale in facies above and below OWC (Figure 3), and constrain its positioning by the transition matrix P (equation 12). However, the present application focuses on identifying the presence of oil and OWC positioning.

There are possible challenges for the inversion because some geophysical rock properties of the same facies are barely correlated between wells, i.e., shale $I_s$ and ${\mathit{Res}}_v$ between WPs and WT; also, some properties of distinct facies show superimposed distributions, i.e., $I_p$ of brine and oil sands for the three wells.

Rock-physics calibration

The rock-physics modeling parameters (see Appendix B and Figure 5) are inverted from the facies-related geophysical and petrophysical properties well logs of WP1 and WP2, using deterministic nonlinear inversion (Tarantola, 2005) and starting from typical values for this environment.

The mean errors in the adjustment are approximately 2.4%–3.8% for WP1 and 3.4%–6.5% for WP2, depending on the facies. Figure 5 shows that most of the transformed elastic properties lie inside the hard data distributions except for WPs shale, due to the few samples and discrepancies in density profiles that lead to a poor fitting of densities, although not of velocities. The impedances $I_p=\mathit{Dens}·V_{\mathrm{p}}$ and $I_s=\mathit{Dens}·V_{\mathrm{s}}$ are products of the actual modeling properties (see Appendix A), only used for plotting.

Figure 5 shows almost perfect fittings between measured and transformed resistivities because the water saturation well log is calculated using the resistivity well log. The contours are unskewed because, in the present approach, ${\mathit{Res}}_v$ is uncorrelated to $S_w$ in the prior covariance matrices because they are properties of distinct kinds. The following methodology evaluates the most relevant aspects of these calibrations in light of the proposed inversion.

Facies from priors

The Bayesian approach is applied to infer the facies model from the synthetic model of geophysical rock properties under different priors as a quality control to support the interpretation of the subsequent geophysical modeling and inversions. Using the averages and covariances of each WP and scanning $\alpha$ allow us to evaluate the correspondence among the priors and the synthetic model, the calibration of the petrophysical modeling, and the linearization.

Only the geophysical rock properties of the synthetic model are used as input for simulating geophysical data in further inversions. Thus, Bayes’ rule for estimating facies conditional on properties m takes the form

where $p(m|\pi )$ is given by equation C-10 and $p(m)$ is a normalization factor as in equation 2.

Addressing the prior probability of facies occurrence and transitions is outside the scope of this analysis; the inference is performed point-wise for each layer (l-index) with a uniform probability of facies occurrence and without setting any correlation between the layers. Thus, substituting equation C-10 in equation 18 produces

where the index i relates to facies ($i=1$ is shale, $i=2$ is brine sand, and $i=3$ is oil sand).

Hence, the MAP facies for each layer is determined by maximizing the facies probability for a given $m_l$:

Figure 6 shows the facies models inferred by the described method from the WPs, from elastic, electric, and joint properties, for 11 $\alpha$ values regularly spaced in the $0\le \alpha \le 1$ range. Note that the IOI top now is rounded to 2800 m and the bottom is 2870.2 m, an interval corresponding to 13 layers of 5.4 m (see Appendix A).

The results show that WP1 priors overestimate the oil column in elastic, retrieve the reference facies in electric, and retrieve the reference oil column in joint plots. The WP2 applications also retrieve the reference facies in electric and the reference oil column in joint plots. In opposition to WP1, the elastic plot in WP2 shows an underestimated oil column. These results support the worth of using resistivity to map fluid variations.

In most cases, the resulting facies show weak dependence on $\alpha$, which validates the proposed adaptation of nonlinear rock-physics modeling to the linear Gaussian approach (equations C-7, C-6, and 10).

Although some of the preceding results present brine just above oil, constraining facies transitions through equation 12 will avoid such a geologic inconsistency in subsequent inversions.

Geophysical linearization

In the “Statistical model” subsection, the geophysical modeling is assumed to be close to linearity surrounding the averages ${\mu }_{\eta }$ to achieve equation 10. Thus, the geophysical responses for the synthetic model of m properties are generated by the nonlinear functions for seismic $g^S(m_S)$ and CSEM $g^{\mathrm{EM}}(\mathit{Res})$ modeling as described in Appendix A and compared with facies-related linear functions $g_{\mathrm{Lin}}^S(m_S)$ and $g_{\mathrm{Lin}}^{\mathrm{EM}}(\mathit{Res})$. In that analysis, the prior means and covariances are extracted from the well logs of WT itself (see Figures 4 and 5). It is limited to geophysical rock properties m ($\alpha =1$ in equation 10) for assessing only linearization effects in geophysical modeling.

As discussed in Appendix A, despite using Zoeppritz for modeling the reflection coefficients of seismic data, the Jacobian uses the Aki-Richards approximation (Aki and Richard, 2002) to seismic data, as follows:

with $m_S$ as the logarithm of the elastic properties and $G_S=S*AD$, where $G_S$ is a function of ${\mu }_{V_{\mathrm{p}},V_{\mathrm{s}}|\pi }$.

Meanwhile, CSEM linear modeling uses the exact analytical derivatives of nonlinear modeling (see Crepaldi et al., 2011 and Appendix A), evaluated at ${\mu }_{\mathrm{Res}|\pi }$ as follows:

Figure 7 compares the geophysical responses of the linear with the nonlinear modeling over the synthetic model for clear and contaminated signals. The difference between the responses appears to be below the noise level for all data sets. However, the 5% difference in the inline electric field at far offsets supports using the nonlinear function even in this case of 1D simulated data.

To compare the linear cases with the proposed nonlinear adaptation in $p(d|\pi )$ of equation 10, Figure 8 shows histograms for acceptance rates (see the “Facies from geophysical data” subsection):

of trial facies models $p(d|{\pi }^{\mathrm{trial}})$ against the reference facies model $p(d|{\pi }^{\mathrm{reference}})$. The histograms represent the cumulative acceptance rates divided into 10 intervals between 0 and 1 for 100 simulated data d with different noise seeds in pictures for seismic and CSEM stand-alone and joint data for 13 model trials.

It is worth pointing out that seismic data uses the Aki-Richards approximation for the simulated data $d_S$ and forward-modeling $g_{\mathrm{Lin}}^S({\mu }_S)$ in the linear case and Zoeppritz for the simulated data and forward modeling in the nonlinear case. In contrast, only the linearly generated $d_{\mathrm{EM}}$ is different from the nonlinear case, whereas the forward modeling is the same for the cases $g_{\mathrm{Lin}}^{\mathrm{EM}}({\mu }_{\mathit{Res}})=g^{\mathrm{EM}}({\mu }_{\mathit{Res}})$ (equation 22).

These histograms for the probability of accepting trials against the reference facies model are proposed as a simple way to assess the insertion of nonlinear modeling on the linear Gaussian approach because conceptually, it should be $r=1$ for the correct trial model and $r=0$ for the others. Although far from sampling the entire domain as the Bayesian inversion intends to do, comparing the posterior likelihood function for linear and nonlinear modeling and the resolution of each geophysical data set in this way is enlightening.

The behavior of the histograms is generally similar for these functions, validating the nonlinear adaptation. Furthermore, the higher sensitivity to contrasts of the nonlinear modeling provides a more accurate acceptance rate for all of the displaced distributions in seismic and CSEM stand-alone and joint data.

These results previously clarify that seismic difficulty in discriminating oil from brine sand and CSEM’s tendency to overestimate the oil column on subsequent inversions are unrelated to adapting to the nonlinear forward functions on the linear Gaussian approach.

Despite these results validating this application, performing similar numerical tests before running the inversion in other scenarios is recommended.

Inversion for facies

All inversions start from a homogeneous model of shale, the facies occurrence probability is uniform ($p({\pi }_1)=\{1/3\}$ in equation 12), and the transition matrix P imposes a 90% chance of choosing the same facies on the layer just below each one, 10% of changing, and zero of placing oil just below a brine layer.

Main application

Figure 9 shows the MAP facies models inverted from contaminated geophysical data with a noise seed for each of 11 realizations using $\alpha$ regularly spaced values that cover the $0\le \alpha \le 1$ range. These realizations are plotted side by side in frames of seismic (SS) and CSEM stand-alone (EM) and joint seismic/CSEM (SEM) inversions. The MAP facies model is the most accepted configuration in 2000 iterations, one order higher than the approximately 200 iterations of the burn-in period.

When comparing the inversions, SS applications show unstable and inaccurate results, contrasting with the EM. because the seismic signal is sensitive to the model contrasts, or differentiation, the CSEM signal is sensitive to the model integration (concept well expressed by the use of transverse resistance in Alvarez et al. [2018], Miotti et al. [2018], and Menezes et al. [2021]), making it less susceptible to uncorrelated noise fluctuations. On the other hand, the integrative aspect of the CSEM signal impacts its resolution, tending to overestimate the oil column by stretching it up and down, also related to WP and $\alpha$.

Fortunately, the SEM results meet the expectations about the complementarity of the geophysical methods when the contribution of the seismic data leads to more stable estimations than EM for the oil column, with an average bias of one extra oil layer below the reference OWC (Figure 9), regardless of the prior settings. That means more reliability for reservoir interpretation and drilling risk assessment.

Lowest resistivity contrast

Multiplying the resistivity profiles of the three used wells by the same anisotropy coefficients ${\theta }_{\pi }$ might be an artificial advantage of inverting simulated data. Hence, similar inversion experiments are performed with the original resistivities to assess the results without such an effect.

In Figure 10, inversion experiments that use the original horizontal resistivities of WPs and WT (Figure 4) present similar behavior to those with high anisotropy coefficients on the oil sand (${\theta }_3=10$) shown in Figure 9, i.e., more stable OWC mapping in SEM than in EM inversions.

Although the higher the oil column resistivity, the higher the CSEM response, applications of vertical resistivity in Figure 9 retrieve worse OWC mapping than in Figure 10. This happens because the measured horizontal electric signal produced by the horizontal electric dipole is dominated by the transversal magnetic mode (Nabighian, 1987), which has a waveguide effect that concentrates the energy in the most resistive layers, shielding the sensitivity to the layers just below the conductive layers, as observed in Figure 8. Aiming to support this argumentation, Figure 11 shows the 13 by 13 elements of sensitivity matrices (Tarantola, 2005) ratio $(G_c^T{G}_c)/(G_r^T{G}_r)$, where $G_c$ and $G_r$ are the CSEM Jacobians of the more conductive and the more resistive models, respectively. Note that the yellow squares (the ratio above one) are only related to the layers below the OWC, whose sensitivity is higher in the more conductive model.

These accurate SEM estimations for the oil column, ten times lower than the typical vertical resistivities in that environment (approximately 70 Ωm) (Zerilli et al., 2018; Menezes et al., 2021), indicate that the proposed methodology might work similarly to the 3D body CSEM response, which also is lower than the 1D modeling used (Crepaldi et al., 2011).

Misguided overburden

Despite the prior and reference models based on well logs and the realistic noise added to the synthetic geophysical data, these controlled inversions have typical advantages over in-field applications: 1D isotropic approximations, exact acquisition parameters, and the same IOI and discretization as the modeled data benefit these geophysical methods. The same background model and anisotropy ratio benefit CSEM, whereas the same wavelets and the convolutional model approximations (see Appendix A) benefit seismic data.

The physics of seismic wave propagation allows the separation of the amplitude-variation-with-angle target signal from the background-related signal by data processing (Castagna and Backus, 1993). In contrast, electromagnetic diffusion is highly influenced by geoelectric background features, with nonlinear effects often inseparable from the target signal. Hence, we must insert the background model into the simulated data and the forward modeling (see Appendix A). This section simulates a challenging scenario for CSEM inversion.

Figure 12 shows a similar experiment. However, the CSEM simulated data are perturbed by a 10 Ωm–20 m layer, inserted 1000 m below the sea floor and 480 m above the IOI that is not introduced in the forward inversion modeling, to simulate a pessimistic scenario in which the CSEM inversion model might bias the interpretation of the reservoir due to the presence of some unidentified resistive rock in the overburden (Zerilli et al., 2018), e.g., carbonates, gas hydrate, salt, igneous and volcanic rocks, and cemented sandstone.

In contrast with the previous results, EM and SEM inversions retrieve the oil column consistently overestimated by four to six extra oil layers below the reference OWC, showing that this anomalous resistive layer strongly impacts the CSEM signal and could introduce errors in its interpretation. However, penalizing the CSEM data fitting in SEM by incorporating a multiplicative factor to increase the corresponding variance elements in $C_d$ (equation 11) allows the seismic signal to improve the OWC estimation (last frame in Figure 12).

Under the Bayesian approach, prior covariances control the weight of each rock property and geophysical data in the facies inversion and could be interpretation parameters, as in the present case. For instance, Gao et al. (2012) propose balancing electromagnetic to seismic data by the ratio between starting model responses misfit.

Such an example suggests that the proposed simultaneous geophysical inversion could mitigate other sources of errors from each geophysical method.

Inversions for rock properties

Figures 13 and 14 show conditional distributions of geophysical and petrophysical rock properties built by sampling deterministic nonlinear inversions performed over the accepted facies models, as described in the “Local inversions” subsection. Figure 13 uses WP1 as a prior and Figure 14 uses WP2 for $\alpha =0.5$.

Figure 15 shows the same inversion experiment with $\alpha =0.5$, although using profiles from the WT itself to gather the prior distributions and calibrate the rock-physics modeling. As expected, these more likely prior inputs provided better estimations than those using WPs.

Figure 16 shows the misfit for each facies iteration on the SEM-WP1 inversion of Figure 13 after the local inversions for geophysical $(d-g[m_k]{)}^T{C}_d^{-1}(d-g[m_k])$ and petrophysical $(m_k-t[r_{k_r}]{)}^T{C}_{m_r}^{-1}(m_k-t[r_{k_r}])$ rock properties. Figure 17 shows the respective geophysical data, sampled, and MAP responses.

The inversions for rock properties are presented as an extra application because it is an easily accessible by-product of the proposed facies inversion. However, approaching that at the same level of detail as the facies inversion would require many more experiments beyond the scope of this work.

DISCUSSION

The run time for a single-loop facies inversion is approximately 7.4 min using 17% capacity of an Intel i9-9900k 3.6 GHz processor with all functions in noncompiled Matlab codes. Including the local inversions for rock properties increases the run time to 56 min under the same conditions.

It is necessary to implement the codes in a compiled language and use parallelization to apply the routine in production for interpreting 2D or 3D bodies, even using the CSEM 1D modeling in common midpoint gathers as Mittet et al. (2008), Buonora et al. (2014), and Corrêa and Régis (2017) do in exploration grids if this approximation also is valid at the reservoir scale. This is because 2D and 3D reservoir models may have hundreds to hundreds of thousands of horizontal cells that must be laterally constrained.

Some other possible advances in this workflow are: (1) optimizing $\alpha$ to condition the inverse matrices in facies and local inversions, (2) postinversion optimization for the OWC, (3) incorporating 3D CSEM and seismic modeling (i.e., finite difference), and (4) including anisotropy as a rock property.

CONCLUSION

The semianalytical Bayesian algorithm successfully estimates models of lithofluid facies conditioned on seismic and CSEM data, including petrophysical and geophysical property distributions, weighted by an introduced parameter that allows us to test the inversions on many prior conditions.

The proposed generalization that associates two prior distributions in the same marginalization integral leads to an expression that prioritizes the most accurate information between them in the posterior covariance matrix. Such an aspect of this original approach also might be attractive for other applications.

The facies inversion is applied to well-based synthetic models with priors from one of two other wells per experiment, simulating scenarios of minimal information compared with what is usually available during the oil field recovery phase. Given that and the added noise, the results show that joint inversion for lithofluid facies can map vertical fluid variations in reservoir scale and robustly deal with noise and geologic uncertainties in prior settings.

In stand-alone applications, the CSEM inversions provided stable oil column estimations inside subsets of prior information, in contrast with the general instability of the seismic inversions. However, the seismic influence on joint inversions stabilized the OWC positioning, providing results weakly dependent on prior settings and data noise, always indicating the presence of oil with an average bias of two oil layers in addition to the six of the true model, with 5.4 m thickness each.

By comparing the inversion results with the incompatibilities between priors and elastic synthetic models previously observed by quality controls, it is found that the remaining deficiency of joint inversion to precisely map the oil column is CSEM resolution limitations and weak elastic response for fluid variations.

Two possible challenging scenarios for CSEM data inversion also are simulated. One with lower resistivity contrasts for fluid saturation provides better OWC positioning because the less resistive oil column weakly shields the sensitivity to the lower layers.

In another scenario, the CSEM simulated data are perturbed by an anomalous resistive layer inserted in the overburden model not included in the inversion background model, leading to the CSEM inversion overestimating the oil column, as expected. However, joining the weight-enhanced seismic data improved the oil column estimation. That represents an example of difficult-to-access compensation for collaborative interpretation or multistep joint inversion, justifying the simultaneous joint inversion and indicating that it can reduce other sources of uncertainty in stand-alone signal interpretations.

Assembling geophysical methods, rock properties, and prior information on the proposed semianalytical facies inversion allows low-cost estimations of conditional distributions for geophysical and petrophysical rock properties. Although applied as a supplementary subroutine without testing many parameterizations, the results demonstrate potential for reservoir quantitative analysis, the quality of which depends on prior inputs.

DATA AND MATERIALS AVAILABILITY

Data associated with this research are confidential and cannot be released.