1. Introduction

Integration of multidisciplinary data represents a fundamental task in geophysics. Nowadays, modern computing technology allows several approaches for combining different types of data, including seismic, electromagnetic, gravity, magnetic and borehole data. Simultaneous Joint Inversion (SJI) represents one of the most sophisticated approaches for optimizing the integration workflow. The final objective is to combine complementary geophysical methods in a single inversion scheme, honoring the total information simultaneously. This can be achieved by minimizing a ‘joint cost function’ including two or more misfit terms, plus one or more regularization terms and, eventually, some type of constraint. For example, in hydrocarbon industry, electromagnetic and/or gravity data are inverted simultaneously with seismic data, producing a consistent model of density, electric and seismic properties. More complex inversion algorithms are used for retrieving 2D or 3D spatial distribution of rock properties like porosity and fluid saturation. Many authors have described different approaches to perform SJI in various fields (Vozoff and Jupp, 1975; Haber and Oldenburg, 1997; De Stefano et al., 2011; Gallardo and Meju, 2003; 2004; 2007; Colombo et al., 2007; Colombo and De Stefano, 2007; Commer and Newman, 2009; Hu et al., 2009; Abubakar et al., 2011, Liang et al., 2011; Moorkamp et al., 2011). Despite the multitude of approaches, there are many common features in all the SJI proposals. Complementary datasets are inverted simultaneously but not independently. In fact, complementary geophysical methods are integrated in a single inversion scheme honoring all data simultaneously. In order to fit all the multidisciplinary data, a “composite misfit functional” (including all the differences between observed and predicted responses in the various geophysical domains) must be minimized. Moreover, the different domains are reciprocally connected using multiple categories of links like empirical relationships, rock-physics models, cross-gradient constraints or model covariance matrix. The unquestionable practical results achieved by the authors in using SJI, are often customized to the particular case history, and it becomes difficult to derive hints for a general SJI stand-alone processing strategy. Inappropriate links between the different types of data, inadequate treatment of the variable resolution of the different data sets, and many other issues, can create artifacts and wrong models. For instance, Conjugate Gradients or Gauss-Newton optimization methods can drive the objective functional to local minima. This is frequently unavoidable, because the huge amount of data to be inverted simultaneously prevents any possibility to explore the model space extensively by global searching methods. The undesirable result is that the final solution is intrinsically non-unique. Moreover, it is often difficult to quantify the uncertainties and the effective reliability of the final multi-parametric models.

Our rationale is that SJI works in the best way if it is run on a starting model that is progressively updated going through complementary integration approaches. These can include single domain inversion, constrained inversion, cooperative modeling and inversion, all of them supervised by a multidisciplinary team. In order to assist this operation mode, we embed our SJI algorithms in a general integration platform, the Quantitative Integration System (QUIS), strongly based on a parallel computing technology. The platform permits an efficient creation of workflows, data and model management, input-output visualization, selection of links between different algorithms and software, interaction between people in collaborative environments. Due to the complexity of all the above discussion points, this paper is divided in two main parts. In PART I we discuss the theoretical aspects of our integrated approach, which is focused on the quantitative estimation of relevant rock physics properties (mainly porosity and fluid saturation). Moreover, we explain how and why our forward and inverse algorithms work better if they are included in a more general integration platform (QUIS). PART II is dedicated to examples by which we intend to clarify our methodology of inversion/integration and its applicability in the real world. We start discussing a 1D synthetic test of SJI. Finally we show a real application of our QUIS platform where integration of multidisciplinary data sets is performed combining cooperative iterative modeling and inversion.

PART I: THEORY

2. SJI of electromagnetic and seismic data

We present an overview of the Simultaneous Joint Inversion problem and of our proposed solution, with particular focus on the combination of electromagnetic and seismic information: the main objective is to obtain consistent multi-parametric models of porosity and fluid saturation in the reservoir rocks by Simultaneous Joint Inversion (SJI) of P-velocity and resistivity distribution.

We focus our attention on Marine Controlled Source Electromagnetic method (CSEM), due to the importance of this prospecting approach, recently introduced in hydrocarbon exploration (Eidesmo et al, 2002). Furthermore, we discuss in detail our SJI algorithm and how it is included in our QUIS platform. In the next two sections we recall just the key elements about electromagnetic and seismic inversion in their specific domains, addressing the reader to our previous works and/or to the extensive literature available about both subjects (Dell’Aversana, 2014).

2.1 The CSEM inversion problem

The CSEM surveying technique investigates the subsurface below the sea bottom by means of low frequency electromagnetic waves (Constable and Srnka, 2007). The dynamics of the transmitted and measured electric field (an analogue equation is valid also for the magnetic field) is described by the damped wave equation:

$\nabla \times \nabla \times \vec E\, \, \, + \, \, \, \mu \, \sigma \frac{{\partial \vec E}}{{\partial t}}\, \, + \, \, \mu \, \varepsilon \frac{{{\partial ^2}\vec E}}{{{\partial ^2}t}} = 0$

(1)

where $\vec E$ is the electric field, $\sigma$ the electric conductivity, $\varepsilon$ the absolute dielectric constant and $\mu$ the magnetic permeability. The second term in Eq. (1) is the loss term, disappearing in free space, while the third term is responsible for wave propagation. For a typical CSEM scenario (frequency smaller than 1 Hz, conductivity between 1 and 10^-6 S/m, dielectric constant between 10^-9 and 10^-11 F/m, $\mu$ constant) the propagation term can be neglected, and Eq. (1) reduces to

$\nabla \times \nabla \times \vec E\, \, \, + \, \, \, \mu \, \sigma \frac{{\partial \vec E}}{{\partial t}}\, \, = 0$

(2)

The EM propagation is now governed only by the electric conductivity $\sigma$ of the medium, or by its reciprocal, the electric resistivity. Oil and gas accumulation can produce resistivity anomalies with respect to the mean properties of the background, and they are visible like perturbations of the electromagnetic response. In summary, the objective of the CSEM technique is the inversion of the subsurface resistivity spatial distribution, starting from the EM responses, and the identification of resistive targets (e.g. hydrocarbon reservoirs).

2.2 The seismic inversion problem

The final goal of seismic inversion is the identification of the elastodynamic properties of a given subsurface domain, starting from the traces collected with the seismic acquisition. Various sets of independent parameters can be considered for the full elastic characterization of the domain. The more common set contains the density $\rho$, and the P and S wave velocities, $\left( {{V_P}, {V_S}, \rho } \right)$. The equations

$\begin{array}{l} {V_P} = \sqrt {\frac{{\lambda + 2\mu }}{\rho }} = \sqrt {\frac{{k + 4\mu /3}}{\rho }} \\ {V_S} = \sqrt {\frac{\mu }{\rho }} \end{array}$

(3)

link the seismic velocities to the bulk modulus $k$ and to the Lamè parameters $\mu$ and $\lambda$, so that new domain parameterizations can be defined using the variable sets $\left( {\lambda , \mu , \rho } \right)$ or $\left( {k, \mu , \rho } \right)$.

The determination of the elastodynamic parameters in a Full Wave Inversion (FWI) analysis is based on the comparison between the measured seismic traces and the corresponding synthetic traces obtained by a full-wave simulation, through the minimization of a given objective function.

2.3 The Simultaneous Joint Inversion problem

The Simultaneous Joint Inversion (SJI) of CSEM and seismic data has the final aim of obtaining a characterization of the subsurfacestructure for both electrical and elastodynamic parameters. It allows reducing the uncertainties which affect the inversion results when running single domain procedures. The key point is the existence of cross-domain relations, that we select at the level of rock properties, where some of them (cross-properties) affect many geophysical domains. As described in the next section, a set of empirical models exist in rock physics, providing the required cross-properties relations.

A specific case of SJI is when electric and magnetic CSEM data with different frequenciesare inverted simultaneously for retrieving a robust resistivity model at variable spatial scale. With this approach, high frequency data are more appropriate to constrain the resistivity of thin layers, like hydrocarbon filled reservoirs. The geometrical properties of the reservoir are commonly fixed using seismic information, whereas the CSEM multi-frequency data are used only for retrieving the resistivity distribution inside the reservoir. Instead, low frequency CSEM data constrain better the resistivity of the geological background. An example of this type of SJI of multi-frequency CSEM data constrained by seismic data will be shortly discussed in the second part of this paper.

2.4 Constitutive equations and rock cross properties

Several models are available in the literature connecting rock properties with geophysical measurements (Schön, 2011). The fundamental remark is that different observations (multidomain geophysical data) converge to a common set of physical parameters (rock properties), making it possible and advantageous the joint optimization of the subsurface imaging problem, with an effective reduction of interpretation ambiguities. For an hydrocarbon exploration scenario, e.g. for the characterization of a porous rock saturated with oil, water and gas, the core parameters are the solids and fluids “mixture” in a geophysical equivalent compound medium, that can be described by the following parameters (suffix “w” indicates water, “o” is oil, “g” is gas, “solid” is the rock matrix):

- bulk modulus of the rock matrix (${k_{solid}}$), and of the filling fluids (${k_w}, {k_o}, {k_g}$);

- shear modulus ${\mu _{solid}}$ of the rock matrix (for the filling fluids μ = 0);

- density of the rock matrix (${\rho _{solid}}$), and of the filling fluids (${\rho _w}, {\rho _o}, {\rho _g}$);

- conductivity of the rock matrix (${\sigma _{solid}}$), and of the filling fluids (${\sigma _w}, {\sigma _o}, {\sigma _g}$);

- porosity $\varphi$;

- water (${S_w}$), oil (${S_o}$), and gas (${S_g}$) saturation, with ${S_w} + {S_o} + {S_g} = 1$.

Seismic geophysical methods recover elastic wave velocities, which in the “core” domain are related to the matrix and fluid bulk modulus, shear modulus, density and porosity. As an example, we show here the Raymer model (Raymer et al., 1980): it combines the velocities of the solid and fluid phases in order to obtain the P-wave velocity of the equivalent medium as

$\begin{array}{l} {V_P} = {(1 - \varphi )^2}{V_{solid}} + \varphi {V_{fluid}}\quad \quad \varphi < 0.37\\ \frac{1}{{\rho V_P^2}} = \frac{\varphi }{{{\rho _{fluid\;\;\;\;}}V_{fluid}^2}} + \frac{{1 - \varphi }}{{{\rho _{solid\;\;\;\;}}V_{solid}^2}}\quad \quad \varphi > 0.47\\ \frac{1}{{V_P^{}}} = \frac{{0.47 - \varphi }}{{0.10}}\frac{1}{{V_{37}^{}}} + \frac{{\varphi - 0.37}}{{0.10}}\frac{1}{{V_{47}^{}}}\quad \quad 0.37 < \varphi < 0.47 \end{array}$

(4)

where ${V_37}$ and ${V_47}$ are the P wave velocities obtained for $\varphi$ equal to 0.37 and 0.47 respectively, and the solid and liquid phase velocities are defined as

$\begin{array}{l} {V_{solid}} = \sqrt {\frac{{{k_{solid\;\;\;}} + 4{\mu _{solid\;\;\;}}/3}}{{{\rho _{solid}}}}} \\ {V_{fluid}} = \sqrt {\frac{{{k_{fluid}}}}{{{\rho _{fluid}}}}} \end{array}$

(5)

with

${k_{fluid}} = {S_w}{k_w}\, + \, {S_o}{k_o}\, + {S_g}{k_g}$

(6)

and

${\rho _{fluid}} = {\rho _w} \cdot {S_w} + {\rho _o} \cdot {S_o} + {\rho _g} \cdot {S_g}$

(7)

Eq. (7) represents the conservation of mass, and it can be extended to the whole compound rock,

$\rho = [(1 - \varphi ) \cdot {\rho _{solid}} + \varphi \cdot ({S_w}{\rho _w} + {S_o}{\rho _o} + {S_g}{\rho _g})]$

(8)

In a similar way we can relate electromagnetic observations, e.g. the conductivity of the equivalent medium as observed with a CSEM acquisition, to the “core” domain properties. The example here are first and second Archie laws. The electric current in saturated porous rocks is almost exclusively due to the conductivity of the saturating fluids, the constituting minerals being semiconductors or insulator in the largest majority of cases. The ratio between the saturating fluid conductivity ${\sigma _f}$ and the rock effective conductivity $\sigma$ represents the formation factor $F$

$F = \frac{{{\sigma _f}}}{\sigma }$

(9)

For a fully saturated rock Archie first law links the formation factor to the rock porosity:

$F = a{\varphi ^{ - m}}$

(10)

where $a$ is an empirical constant with value close to 1 and the exponent $m$ (cementification exponent) is between 1.3 and 2.5 for sedimentary rocks, and close to 2 for sandstones. It must be noted that the first Archie law is only valid for clean rocks: in presence of shaly sands other empirical laws must be used, as detailed below.

The Archie second law is a relation between the effective conductivity of a partially saturated rock (e.g., a rock whose saturating fluids are water and oil in comparable fractions) with the water saturation factor ${S_w}$ and conductivity ${\sigma _w}$,

$\sigma = \frac{{S_w^n}}{F}{\sigma _w}$

(11)

$n$ is the saturation exponent, whose value depends on the saturation conditions (saturating fluids and fractions): rocks for which the fraction of oil is larger than the water one show a larger $n$ value with respect to rocks whose saturating fluid is prevalently water.

In general, constitutive equations are empirical/physical functions g, f and h relating the equivalent medium P-velocity ${V_P}$, electrical conductivity $\sigma$, and density $\rho$ to the porosity $\varphi$ and the vector ${\bf{s}} = {\left[{{S_w}, {S_o}, {S_g}} \right]^T}$ of fluidsaturations,

$\begin{array}{l} {V_p} = {\rm{g}}(\varphi , {\bf{s}}), \\ \sigma = {\rm{f}}(\varphi , {\bf{s}}), \\ \rho = {\rm{h}}(\varphi , {\bf{s}}). \end{array}$

(12)

In this way we are pointing out that porosity and fluid saturations are common parameters of different geophysical domains (cross-properties), and this permits to link the domains with cross-relation of the type

$\sigma = {\rm{f}}\, {\rm{[}}{{\rm{g}}^{ - {\rm{1}}}}({V_p})]$

(13)

The crucial step of selecting the appropriate set of constitutive equation is driven by geological information, and/or by experimental analysis of real data collected in situ.

Furthermore, cross-property relations are useful when some rock properties are obtained more easily than other properties. For example, when seismic velocity can be estimated more accurately than conductivity, the latter can be derived from the cross property relation (13).

Following this approach, from a set of constitutive equations, we can identify a set of cross-properties, and a set of cross-relations.

3. General formulation of Simultaneous Joint Inversion

Our inversion algorithm requires in input the grid data of P-velocity, density and conductivity derived from a different single domain inversion or a priori information. The first step of the inversion process is to map the available domain on a common domain and a common grid. First of all the inversion domain is defined as the intersection of the single domains and the data are mapped on the finer available grid in order to have the geophysical properties in the same points.

The Joint inversion problem is solved by minimizing a Joint Objective Functional (JOF). This can be schematically represented by the sum of the following terms:

· Joint misfit functional,

· Regularization terms,

· Constraints and relationships.

We apply the stochastic approach of Tarantola (2005), where the misfit functional for a given inverse problem assigned in a certain geophysical domain is defined as

${\Phi _{mis}}({\bf{m}}) = {\left[{g({\bf{m}}) - {{\bf{d}}_{obs}}} \right]^T}{\bf{C}}_D^{ - 1}\left[{g({\bf{m}}) - {{\bf{d}}_{obs}}} \right]$

(14)

where m and d_obs are, respectively, the multi-parametric model vector and the data vector. For instance, in case of joint petro-physicalinversion, the multi-parametric model vector can include porosity and fluid saturations, whereas the multi-parametric data vector can include velocity, density and conductivity, g(m) is the simulated response calculated using the forward operator g(nonlinear in the most general case), and ${\bf{C}}_D^{ - 1}$ is the inverse of the data covariance matrix.

For each domain, a specific forward operator is applied to the correspondent model vector, producing a theoretical response to be compared with the observations. The joint misfit function will be simply the sum or a linear combination of the individual misfit functions.

Different types of regularization can be included in the objective function in order to stabilize the inversion results (Zhdanov, 2009). A frequently used regularization approach is based on smoothing operators (Vozoff and Jupp, 1975), like the following (Moorkamp et al., 2011):

${\Phi _{reg}}({\bf{m}}) = \sum\limits_i {{\alpha _i}{{({\bf{m}} - {{\bf{m}}_0})}^T}W_i^T{\bf{C}}_M^{ - 1}{W_i}({\bf{m}} - {{\bf{m}}_0}) + \beta } {({\bf{m}} - {{\bf{m}}_0})^T}{\bf{C}}_M^{ - 1}({\bf{m}} - {{\bf{m}}_0})$

(15)

In equation (15) $i = \left\{ {x, y, z} \right\}$ represents the three spatial directions and ${\alpha _i}$ is a weighting factor; ${{\bf{m}}_0}$ is a reference model; ${\bf{C}}_M^{ - 1}$ is the inverse of the model covariance matrix; and the matrix ${W_i}$ represents an approximation of the spatial (first or second) derivative of the model parameters.

Each single inversion domain point can be decoupled during the inversion process in order to obtain the corresponding optimal petrophysical properties (like porosity $\varphi$ and saturations ${S_w}$, ${S_o}$ and ${S_g}$).

A set of constitutive relations is chosen: potentially, for each point, it is possible to set different constitutive relations and different model parameters,

$\begin{array}{l} {V_p}\, \, = \, \, {V_p}\, \, \left( {\varphi \, , \, {S_\omega }\, , \, {S_o}\, , \, {S_g}} \right)\\ \sigma \, \, = \, \, \sigma \, \, \left( {\varphi \, , \, {S_\omega }\, , \, {S_o}\, , \, {S_g}} \right)\\ \rho \, \, = \, \, \rho \, \, \left( {\varphi \, , \, {S_\omega }\, , \, {S_o}\, , \, {S_g}} \right) \end{array}$

(16)

The petrophysical properties are obtained following a probabilistic inversion process based on the Bayes theorem, allowing the introduction of a priori information obtained independently into the inversion process. The results of the adopted Bayesian inversion procedure will be the optimal porosity and fluids saturation map of the subsurface and also a measurement of the solution uncertainty. It allows also to take into account the model and data uncertainty, and to introduce the a priori information into the inversion process. In the least-squares formulation of the inverse problem, the (twice) objective function to minimize is given by:

$\begin{array}{l} 2S({\bf{m}}) = ||g({\bf{m}}) - {{\bf{d}}_{obs}}||_D^2 + ||{\bf{m}} - {{\bf{m}}_{prior}}||_M^2 = \\ = {\left( {g({\bf{m}}) - {{\bf{d}}_{obs}}} \right)^T}C_D^{ - 1}\left( {g({\bf{m}}) - {{\bf{d}}_{obs}}} \right) + {\left( {{\bf{m}} - {{\bf{m}}_{prior}}} \right)^T}C_M^{ - 1}\left( {{\bf{m}} - {{\bf{m}}_{prior}}} \right) \end{array}$

(17)

It is possible to verify that the a-posterioriprobability density is approximately Gaussian, with a centre given by

${{\bf{m}}_{est}} \cong {{\bf{m}}_{prior}} + {({{\bf{J}}^T}{\bf{C}}_D^{ - 1}{\bf{J}} + {\bf{C}}_M^{ - 1})^{ - 1}}{{\bf{J}}^T}{\bf{C}}_D^{ - 1}({{\bf{d}}_{obs}} - g({{\bf{m}}_{prior}}))$

(18)

Here J indicates the jacobian operator.

The corresponding a-posteriori covariance operator, which describes the uncertainty of the estimated model, is given by

${\bf{C}}_M^{post} = {\left( {{{\bf{J}}^T}{\bf{C}}_D^{ - 1}{\bf{J}} + {\bf{C}}_M^{ - 1}} \right)^{ - 1}} = {{\bf{C}}_M} - {{\bf{C}}_M}{{\bf{G}}^T}{\left( {{\bf{J}}{{\bf{C}}_M}{{\bf{J}}^T} + {{\bf{C}}_D}} \right)^{ - 1}}{\bf{J}}{{\bf{C}}_M}$

(19)

The typical numerical strategy for these problems is to use some iterative algorithm to obtain the maximum likelihood estimate of the model vector. This means that we start from an initial model (starting guess) and then we update it iteratively. Using, for instance, a quasi-Newtonian optimization approach, the iterative formula is

${{\bf{m}}_{n + 1}} = {{\bf{m}}_n} - {\mu _n}{({\bf{J}}_n^T{\bf{C}}_D^{ - 1}{{\bf{J}}_n} + {\bf{C}}_M^{ - 1})^{ - 1}}[{\bf{J}}_n^T{\bf{C}}_D^{ - 1}({{\bf{d}}_n} - {{\bf{d}}_{obs}}) + {\bf{C}}_M^{ - 1}({\bf{m}}_n^{} - {{\bf{m}}_{prior}})]$

(20)

4 Embedding the SJI algorithm in the general integration platform

Owing to the intrinsic limitations of every individual integration approach, including Joint inversion algorithms, we combined different integration techniques (like single-domain, constrained, cooperative, model-driven and joint-inversion techniques) in an optimized platform. This is called QUIS, Quantitative IntegrationSystem, and it represents a sort of “second-order integration”, i.e. an optimized combination of integration methods. In such a way, the problem of integrating different types of geophysical data becomes a question of building an efficient workflow, data and model management, input-output flow, links between different algorithms and software, and interaction between people in collaborative environments. Here we introduce just the crucial aspects of the QUIS concept. Additional details can be found in previous published works (Dell’Aversana, 2014). Figure 1 shows a conceptual scheme of the QUIS.

This platform consists of linked software that can be properly parameterized with the support of a graphic user interface (G.U.I.). The first set (Figure 1.left) includes original proprietary codes. These are dedicated to solve the forward and the inverse problem in three geophysical domains: seismic, electromagnetic and gravity. Modeling and inversion can be performed in separate domains as well as simultaneously, using one or more integration approaches (constrained, cooperative, joint inversion). In fact, all three domains are linked using a library of rock physics relationships (indicated in figure as “QUIS core based on rock physics”). This “library” contains all the most important models in literature (Schön, 2011), such as the models of Raymer, Archie, CRIM, Hashin-Shtrikman, Gassmann, Biot-Geertsma-Smit, Eberhart-Philips, Eshelby, Wood, Glover, Hermance, Koesoemadinata and McMechan’s cross models, Faust cross model. Moreover it is continuously updated with empirical relationships retrieved from composite well logs.

For instance the user can decide to combine the Raymer model with the CRIM model for running a joint inversion of seismic velocity and resistivity data with the objective to retrieve a porosity/saturation model (see example in the next section).

This can be done in 1D (example: joint inversion of composite well logs), in 2D and in 3D (example: joint inversion of surface geophysical data). All these codes can cooperate with pre-existing software platforms already available. This can happen by exchanging input data and output models in the proper formats, in order to include also “standard” processing operations in the whole integration work-flow. For instance, pre-stack depth migration of seismic data can be performed using a commercial software and the output can be imported in the QUIS platform for triggering a loop of cooperative inversion between seismic and electromagnetic data. Then, the results can be used as a starting model for triggering a simultaneous joint inversion of seismic and electromagnetic data. The final output will be a multi-dimensional porosity/saturation model (see the first example discussed in the PART II of this article).

The fundamental characteristic of this systemic approach is that joint inversion represents a crucial part of the integration work flow, but it doesn’t work as a stand-alone tool. It is well linked with other integration approaches and it is based on criteria of geophysical/geological consistency and on rock-physics relationships.

In the following sections we discuss a synthetic example of simultaneous joint inversion with the purpose of showing what type of results this approach can produce. Moreover, we discuss a case history where our QUIS platform is applied to a real exploration problem.

PART II: examples

In this second part we discuss two examples: the first is a synthetic test aimed at clarifying our joint inversion approach; the second is a complex case history aimed at describing how we combine multidisciplinary data sets in a real exploration project using our QUIS approach.

5 Synthetic test: multi-resolution and multi-domain joint inversion

Our joint inversion framework is validated though a synthetic example that makes use of the elastic full waveform inversion together with the CSEM inversion technique.

The scenario is a reservoir characterization, exploiting prior information coming from different datasets and geological interpretations. The positions of main seismic interfaces are well known, but we have only a limited knowledge of geological formations and lithology from well logs acquired in the prospect area. Moreover, we make the hypothesis to have limited information about petrophysical properties of the reservoir, such as oil saturation and porosity.

The key idea is to combine the information coming from CSEM and seismic. Electromagnetic fields bring information on the subsurface resistivity distribution, also related to water/oil saturation and effective porosity. On the other hand, seismic waveforms are influenced by acoustic impedances, more sensitive to gas-liquid contrasts, rather than to water-oil interfaces.

As previously explained, electromagnetic and elastic/acoustic data can be linked by means of a petrophysical set of equations in order to jointly solve liquid-gas and water-oil ambiguities.

The conceptual scheme of the example can be divided in two main steps (Figure 2): the forward modeling and the inverse problem solution. These two steps can be split in sub-steps.

In the forward modeling phase, a petrophysical “true” test model is built combining information coming from different datasets and geophysical domains. Elastic and CSEM simulators process the petrophysical model yielding seismic and EM data. The joint inference framework inverts synthetic “real” data in order to estimate the inverted model to be compared with the “true” one.

a) Forward modeling

The forward simulation phase produces geophysical electromagnetic and seismic data for a given geological model, described by its petrophysical variables.

The geometrical model is intentionally simple, in order to focus the attention on the geophysical problem, neglecting all dimensional and computational issues. It consists of a sea bottom sequence of plane-parallel layers, representing a 100 m thick sand reservoir embedded in shale formations (Figure 3).

The 100 m thick reservoir is embedded between shale layers and it is composed of nine plane parallel layers, with 11.1 m of average thickness h, 0.24 of average porosity and 0.50 of oil saturation.

The sand reservoir is divided into nine thin layers of different thickness, in order to mimic a transitional porosity and fluid saturation; the fluid saturation varies from water (bottom) to oil (upper layers). In this way the seismic inversion becomes intentionally insensitive to the fluid contents.

We also make the hypothesis that the electrical properties of sands are governed by the Archie's law, the seismic properties by the Raymer equation, and the density by a weighted average model: this set of cross-domain equations are handled by the QUIS core library.

All physical parameters involved in the elastic and in the electrical petrophysical core are listed in Table 1. The constitutive model is expressed in terms of solid, fluid contents and empirical coefficients.

The parameters for the seismic modeling are compressional waves velocity V_P, shear waves velocity V_S and density , as a combination of fluid and solid matrix,

${V_P} = \sqrt {\frac{{k + 4/3\mu }}{\rho }}$

(21)

${V_S} = \sqrt {\frac{\mu }{\rho }}$

(22)

$\rho = \phi {\rho _{fluid}} + \left( {1 - \phi } \right){\rho _{solid}}$

(23)

where $\phi$ is the effective porosity, $\rho$_fluid and $\rho$_solidrespectively fluid and solid densities.

The pressure velocities in the various specimens can be approximated by means of the Raymer equations:

${V_{solid}} = \sqrt {\frac{{{k_{quartz}} + 4/3{\mu _{quartz}}}}{{{\rho _{quartz}}}}}$

(24)

${V_{fluid}} = \left[{{S_o}{V_{oil}} + \left( {1 - {S_o}} \right){V_{water}}} \right] = \left[{{S_o}\sqrt {\frac{{{k_{oil}}}}{{{\rho _{oil}}}}} + \left( {1 - {S_o}} \right)\sqrt {\frac{{{k_{water}}}}{{{\rho _{water}}}}} } \right]$

(25)

Assuming also, for sake of simplicity, a fixed relation between V_P and V_S , the Raymer equation gives the velocity of the compound medium as

$\left\{ {\begin{array}{*{20}{c}} {{V_P} \simeq {{\left( {1 - \phi } \right)}^2}{V_{solid}} + \phi {V_{fluid}} = {{\left( {1 - \phi } \right)}^2}\sqrt {\frac{{{k_{quartz}} + 4/3{\mu _{quartz}}}}{{{\rho _{quartz}}}}} + \phi \left[{{S_o}\sqrt {\frac{{{k_{oil}}}}{{{\rho _{oil}}}}} + \left( {1 - {S_o}} \right)\sqrt {\frac{{{k_{water}}}}{{{\rho _{water}}}}} } \right]}\\ {{V_S} = {V_P}/\sqrt 3 } \end{array}} \right.$

(26)

Analogous considerations can be applied to density, expressing the fluid contribution independently:

${\rho _{fluid}} = {S_o}{\rho _{oil}} + \left( {1 - {S_o}} \right){\rho _{water}}$

(27)

and writing the compound medium density in the following way:

$\rho = \phi {S_o}{\rho _{oil}} + \left( {1 - {S_o}} \right){\rho _{water}} + \left( {1 - \phi } \right){\rho _{solid}}$

(28)

For a two-phase saturating fluid, Archie’s law can be rewritten with respect to oil saturation:

$resistivity = \frac{a}{{{\sigma _{water}}}}{\phi ^{ - m}}{\left( {1 - {S_o}} \right)^{ - n}}$

(29)

Equations (24), (25) and (26) represent the petrophysical parameterization of seismic quantities: they transform the petrophysical space into the geophysical space and viceversa (Figure 4), providing the possibility to describe the whole model in terms of porosity, fluid saturation and constitutive parameters.

The synthetic elastic response is calculated by means of the reflectivity method (Thomson, 1950), with the parameters listed in Table 2.

The CSEM response, namely the electrical and magnetic complex fields E_x and H_y, is computed with a proprietary 2.5D finite element simulator (base on Jin, 2002), with the parameters in Table 3.

b) Inversion

The inversion procedure follows the well-known Bayesian framework, where prior and posterior states of information for data and model parameters are described by Gaussian probability density functions.

The posterior state of information, namely the solution, is obtained after the convergence of a quasi-Newton iterative scheme, in which the a-priori acts simultaneously as a weak boundary and as regularizer (equations 18 and 19).

The key idea is to overcome the limit imposed by the difference of resolutions between seismic methods and CSEM, implementing a novel cooperative and joint scheme of inversion.

We design a workflow of three sequential inversions.

First, we perform a high-resolution single domain inversion using seismic data only in the petrophysical space porosity-saturation. Then, this solution is downscaled and used as a-priori in a joint seismic/CSEM inversion. Finally the new solution is interpolated and used as prior information for a refined seismic inversion.

(i) High resolution seismic full waveform inversion in the petrophysical space

In this step, we discretize the reservoir in ten regularly sampled layers, whereas the layer thicknesses of the “true” model are not regular; this introduce a done-on-purpose systematic noise in the problem, due to the misaligned grids.

The seismic inversion problem is solved in the petrophysical space ($\phi$, S_o) exploiting the parameterization formalized in Equation (26). Each layer of the 1D model is described by a couple porosity-saturation, properly indexed in the vector of model parameters m:

$m = {\left[{{m_\phi };{m_{{S_o}}}} \right]^T}$

(30)

$\left\{ {\begin{array}{*{20}{c}} {{m_\phi } = \left[{{\phi _1};{\phi _2};{\phi _3}; \ldots ;{\phi _M}} \right]}\\ {{m_{{S_o}}} = \left[{{S_o}_1;{S_o}_2;{S_o}_3; \ldots ;{S_o}_M} \right]} \end{array}} \right.$

(31)

where [.]^T is the transpose operator, and M the number of layers.

The data is represented in vector way, collecting seismic samples s of all seismic shot gathers into a vector u with length h

$u = {\left[{{s_1}{s_2}{s_3} \ldots {s_h}} \right]^T}$

(32)

The Jacobian matrix of the partial derivatives of seismic data with respect to the petrophysical parameters is

${G_u} = \left[{\frac{{\partial u}}{{\partial m}}} \right]$

(33)

that can be trivially calculated by a numerical centred finite differences scheme, with respect to perturbations of m.

The chosen prior state of Gaussian information is the following:

$\left\{ {\begin{array}{*{20}{c}} {{m_{prior}} = {{\left[{E\left[\phi \right];E\left[{{S_o}} \right]} \right]}^T}}\\ {{C_M} = diag\left\{ {Var\left( {{m_{prior}}} \right) = {{\left[{Var\left( \phi \right);Var\left( {{S_o}} \right)} \right]}^T}} \right\}} \end{array}} \right.$

(34)

where E[.] is the arithmetic mean and Var(.) the variance operators, both applied to the forward model petrophysical parameters ($E\left[\phi \right] = 0.25, E\left[{{S_o}} \right] = 0.50, Var\left( \phi \right) = 0.031$ and $Var\left( {{S_o}} \right) = 0.15$).

Finally, the iterative Gauss-Newton scheme becomes

${m_{n + 1}} = {m_n} - {\left[{\left( {{G_u}} \right)_n^TC_D^{ - 1}\left( {{G_u}} \right)_n^T + C_M^{ - 1}} \right]^{ - 1}}\left[{\left( {{G_u}} \right)_n^TC_D^{ - 1}\left( {{u_n} - {u_{obs}}} \right) + C_M^{ - 1}\left( {{m_n} - {m_{prior}}} \right)} \right]$

(35)

At the end of the iterative scheme, the posterior state of information is described by its central value and covariance matrix

$C_u^{post} = {\left( {G_u^TC_D^{ - 1}{G_u} + C_M^{ - 1}} \right)^{ - 1}}$

(36)

According to the physics of this problem, we expect that seismic full waveform inversion will be highly sensitive to porosity, but not to saturation due to the weak velocity contrast between oil and water.

The results of petrophysical full waveform inversion confirm the expectations (Figure 5 and Figure 6): in spite of the impressive fitting, sand porosity is resolved while oil saturation proves to have no influence in the seismic data.

(ii) Low resolution joint inversion in the petrophysical space

The solution estimated in step (i) together with its posterior covariance matrix, are used as the prior state of information for the joint seismic/CSEM inversion.

In this step, the inversion grid is downscaled up to 25 m of layer thickness, in order to honor the CSEM resolution at 1 Hz. The posterior probability density function obtained in step (i) refers to a system with 10 layers with thickness of 10m and 20 degrees of freedom; in order to use that information we calculate a down-sampled version with 4 layers and 8 degrees of freedom,

$\left\{ {\begin{array}{*{20}{c}} {{m_{prior}} = {{\bar m}_u}}\\ {{C_M} = diag\left\{ {C_u^{post}} \right\}} \end{array}} \right.$

(37)

Because seismic traces and electromagnetic amplitude/phase recordings are inhomogeneous, data have to be regularized in order to obtain a composite dataset to invert simultaneously.

Inline electric E_x and crossline magnetic H_y field data are pre-regularized using phases in radians, $\phi$_E, $\phi$_H, and the logarithm of amplitudes normalized to the dipole-moment p. In order to obtain the joint regularized data vector, seismic traces and pre-regularized EM data are equalized by means of a multiplier w based on data dynamics.

Pre-regularized EM data are gathered into the vector v with length k

$v = {\left[{{c_1};{c_2};{c_3}; \ldots {c_k}} \right]^T}$

(38)

c is the pair of vectors e and h

$c = {\left[{e;h} \right]^T}$

(39)

Each element of the vectors e and h contains the following values

$\left\{ {\begin{array}{*{20}{c}} {e = {{\left[{\log \left( {\left| {{E_x}} \right|/p} \right);{\varphi _E}} \right]}^T}}\\ {h = {{\left[{\log \left( {\left| {{H_y}} \right|/p} \right);{\varphi _H}} \right]}^T}} \end{array}} \right.$

(40)

The weight coefficient w is

$w = \frac{{MAX\left( v \right) - MIN\left( v \right)}}{{MAX\left( u \right) - MIN\left( u \right)}}$

(41)

with MAX(.) and MIN(.), respectively the data vector maximum and minimum operators.

The composite data vector d becomes

$d = {\left[{wu;v} \right]^T}$

(42)

The electromagnetic Jacobian matrix G_v takes the form

${G_v} = \left[{\frac{{\partial v}}{{\partial m}}} \right]$

(43)

Similarly to the petrophysical seismic inversion, G_v can be calculated perturbing m according to a numerical centered finite differences scheme. Therefore, the joint Jacobian matrix J becomes

$J = \left[{\begin{array}{*{20}{c}} {w{G_u}}\\ {{G_v}} \end{array}} \right]$

(44)

The Gauss-Newton scheme takes the following joint form

${m_{n + 1}} = {m_n} - {\left( {{J^T}C_D^{ - 1}J + C_M^{ - 1}} \right)^{ - 1}}\left[{{J^T}C_D^{ - 1}\left( {{d_n} - {d_{obs}}} \right) + C_M^{ - 1}\left( {{m_n} - {m_{prior}}} \right)} \right]$

(45)

At the end of the iterative scheme, the posterior state of information is described by its central value ${\bar m_j}$ and covariance operator:

$C_j^{post} = {\left( {{J^T}C_D^{ - 1}J + C_M^{ - 1}} \right)^{ - 1}}$

(46)

The fitting of the resistivity model is very good (Figure 7 to Figure 9), considering also the reduced number of layers: the joint utilization of CSEM and seismic data has resolved the saturation ambiguity.

Data are ordered in a regularized super-vector sorted by Ex amplitude, Ex phase, Hy amplitude and Hy phase for each common receiver.

(iii) High resolution seismic inversion refinement in the petrophysical space

In this last step, the same mathematical formulation and parameterization of (i) are used, except for the a priori state of information, which comes from the joint inversion. On the contrary of the previous step, we perform now an upscaling from the 4 layer case to the 10 layer case, by means of linear interpolation:

$\left\{ {\begin{array}{*{20}{c}} {{m_{prior}} = {{\bar m}_j}}\\ {{C_M} = diag\left\{ {C_j^{post}} \right\}} \end{array}} \right.$

(47)

The results of this second seismic inversion show an impressive data fit (Figure 10 and Figure 11): oil saturation has been constrained by the output of the previous joint inversion, and porosity has been updated thanks to a better starting model.

Table 4 is the evolution of the confidence in the inversion result, expressed as the standard deviation of the solution, from the a posteriori model covariance matrix (that becomes the a priori for the following step). Interesting to notice that oil saturation is hardly observable in the seismic data, as the confidence in the “elastic” inverted oil saturation model has a very small improvement. Every step is adding information to the model at the appropriate spatial scale and to the observable parameters, keeping the non observable ones close to the a priori setting. Nevertheless the full sequential flowchart is able to reduce significantly the uncertainty of the final solution.

6 A case history of integration using QUIS

In this section we describe an example of cooperative modelling and inversion of multi-domain geophysical data based on our QUIS approach. In this specific case we combined 3D seismic, electromagnetic, gravity and well data for mapping the reservoirs and for mitigating the exploration risk in a very complex geological setting. This is characterised by narrow and elongated fault compartments with thin stacked reservoir sandstones. The hydrocarbon field has been explored by extensive 2D and 3D seismic surveys, by gravitysurveys and by several wells penetrating hydrocarbon bearing rocks of Upper, Middle and Lower Triassic ages. Permian carbonates have been drilled in the South fault compartment. We recordedmagnetic and electric fields using the marine CSEM method. We acquired a total of 350 km of CSEM towing lines using two central frequencies of 0.5 Hz and 0.125 Hz. All the different types of data (seismic, electromagnetic, gravity and well data) where integrated using our QUIS platform. Figure 12 shows a synoptic view of the complex workflow that we applied.

The workflow starts from the top-right corner of the figure. It continues following the arrows labelled with consecutive numbers. The first step is seismic PSDM, whereas the final step is the transformation of the resistivity models into a probability map of hydrocarbon distribution (step 6).

Due to the complexity of this workflow, we explain the details in the following sub sections.

a) Seismic data interpretation

We started with 3D Kirchhoff pre stack depth migration (PSDM) of seismic data (top right panel in Figure 12). For this step we used a software package external to our QUIS. In fact, our integration platform is optimised in order to combine proprietary algorithms with commercial software packages.

Figure 13 is an example of PSDM section running from NW to SE. The figure shows the interpretation of the top of two reservoirs (in this section indicated as “reservoir 1” and “reservoir 2” or, in the following, as “upper reservoir” and “lower reservoir”, respectively) and of the top of the carbonate platform. A big fault appears towards the SE edge of the section. The reservoirs are both complex and discontinuous, due to the presence of several fault systems, clearly visible in the figure. One of the main difficulties was to understand the lateral extension of the hydrocarbon filled reservoirs. As a general rule, seismic data provide indications about the depth and the geometry of geological layers, but not about the fluid distribution. Instead, CSEM data are sensitive to electric resistivity that, in turns, depends on hydrocarbon saturation. Moreover gravity data allow constraining the interpretation of the main geological units and of the main fault systems. For that reason, we used seismic, CSEM and gravity data as complementary methods for producing a probabilistic hydrocarbon map (as explained in the following).

The direction of this section is showed on the small map. This represents the top of the carbonate platform (in depth) underlying the hydrocarbon reservoirs.

b) Gravity data analysis

After obtaining a seismic cube migrated in depth and after seismic interpretation, we performed the following steps:

■ Bouguer anomaly analysis and iterative forward gravity modelling (step 2 in Figure 12);

■ CSEM data analysis and iterative CSEM modelling (step 3 in Figure 12).

We started with a simple (but useful) co-rendering of CSEM, gravity and seismic information. Figure 14 shows a map of the normalized electric field (for a frequency of 0.5 Hz and an offset of 6.5 km) observed at sea floor superimposed on a high pass filter of the Bouguer anomaly.

Figure 15 shows a detail of the first vertical derivative (FVD) of the Bouguer map (in colours). FVD represents a spatial high pass filter of the Bouguer anomaly. In the figure, it is co-rendered with the Carbonate depth map interpreted from seismic data (contour lines). Figure 16 is an example of FVD profile extracted along the seismic section of Figure 13.

Synthetic gravity data was computed iteratively from a starting model derived from seismic data interpreted in depth (see first step based on PSDM). The initial density model was obtained going through the rock physics module of the QUIS. In fact, we applied different types of empirical relationships in order to retrieve a reasonable density model from the PSDM velocity model.

Gravity information provides the general trend of the main geological structures. This information is useful; however, in order to map the hydrocarbon distribution, we need to analyse the CSEM data. In fact, the main problem remains to map the lateral extension of the oil filled stacked reservoirs. Resistivity is a good discriminator between hydrocarbons (high resistivity) and brine (low resistivity). In order to map the boundaries of high-resistivity anomalies we used a CSEM attribute introduced in previous works. It is given by the normalised difference between out-towing and in-towing CSEM data (Dell’Aversana and Zanoletti, 2010). This is called “CSEM asymmetry attribute”. We demonstrated that it works as a powerful indicator of the edge of resistive bodies. Figure 17 shows an example of this attribute plotted (for two different frequencies) along the seismic section of Figure 13. We interpreted the pick of this attribute as the lateral edges of the hydrocarbon filled reservoir.

c) SJI of multi-frequency CSEM data constrained by seismic data

All the above steps (from “Start” to “step 3” in Figure 12) allowed us preparing a set of robust 2D starting models for the subsequent CSEM inversion work. In order to refine our resistivity models we jointly inverted simultaneously electric and magnetic multi-frequency CSEM data. The inversion was focused on the resistivity parameter, being the geometry of the model already well defined by seismic data. This step (“step 5” in Figure 12) produced a set of 2D resistivity-layered models for the reservoir and for the background. The final models have been updated taking into account the information derived by CSEM asymmetry attribute analysis. This point can be properly understood looking at Figure 18 that combines the previous Figure 17 with the inverted resistivity model.

In the figure we can see that the hydrocarbon filled reservoirs are not laterally continuous. The two gaps, as well as the lateral edges of the reservoirs have been interpreted with the support of the CSEM asymmetry attribute. Instead, the resistivity values have been retrieved by joint inversion of multi-frequency CSEM data (electric and magnetic fields).

We remark that, at this point of the discussion, we are talking about resistors and not about hydrocarbon distribution. Interpreting the nature of resistive anomalies is possible only by combining different sources of information, such as CSEM, gravity and seismic data.

d) From resistivity to hydrocarbon probability maps

The CSEM models provide information about resistivity distribution, but not directly about fluid distribution. Gravity/density information (with the support of seismic data) can be used for reducing the interpretation uncertainty of CSEM anomalies (Dell’Aversana et al., 2011; 2012). The key concept is that the degree of spatial correlation between gravity and CSEM anomalies can be used as a robust indication for interpreting the geological nature of both gravity and electromagnetic responses. For instance, in correspondence of the boundaries of a carbonate formation characterized by high resistivity and high density, we expect to observe high value of both electric and Bouguer anomalies. In that case we observe correlated CSEM and gravity responses at surface. On the other side, the gravity and the CSEM methods have different spatial resolution. In case of relatively small changes at reservoir scale, for instance caused by variations in fluids, we expect to observe only CSEM effects at surface. The gravity method has too low resolution for being sensitive to such a small fluid change. Consequently, in that case, we can see uncorrelated CSEM and gravity responses. Therefore, the level of correlation (or un-correlation) between CSEM and gravity responses can provide useful indications for deciding if an electromagnetic anomaly is more probably caused by presence of hydrocarbons or by resistive geological formations.

Following this criterion we tried to delineate a map of hydrocarbon distribution for each reservoir layer. First we identified the boundaries of the main resistive layers by applying the integrated approach explained above to all the CSEM lines. Where the resistivity anomalies were not explainable in terms of lateral geological variations, they were interpreted in terms of hydrocarbon effects. This result allowed us to define a probabilistic map of hydrocarbon distribution in the exploration block (“step 6” in Figure 12). Figure 19 shows part of the final results. The small map shows the probabilistic distribution map (P) of hydrocarbon in one of the reservoirs (upper reservoir; red means P = 1, blue means P = 0). This map is superimposed on the map of the top carbonates derived from 3D seismic data (large map in background). The parameter P is the probability that the CSEM resistivity anomalies are explained by hydrocarbon filled reservoir (and not by carbonate or other resistive geological units). P is calculated using a Bayesian approach, estimating the posterior probability for a hydrocarbon scenario, in presence of uncorrelated CSEM and gravity anomalies.

In summary, the full integration of all the available geophysical data using the QUIS platform helped to improve reservoir delineation and provided useful indications for mitigating the exploration risk in the area.

7 Conclusions

Geophysical exploration builds images of some physical parameters of the subsurface by generating and recording proper signals. Seismic uses elastic waves, CSEM uses electromagnetic waves, and gravimetric method uses gravity fields. Every domain has intrinsic limitations due to limited resolution, sensitivity and observability. The current trend is to integrate data coming from different domains, exploiting their complementarity, thanks to cross-domain relations between rock parameters. In this framework, many authors are proposing different approaches to perform a Simultaneous Joint Inversion (SJI): the results show big improvements in the rejection of ambiguities and in the increase of the interpretation reliability. Nevertheless, there is not a “standard” flowchart to perform a multidomain multiresolution geophysical data SJI. Thus, we propose to follow a “systemic approach”, and to merge all the software and hardware tools in an integration platform, a Quantitative Integration System (QUIS). The objective of the QUIS is to make available to a multidisciplinary experts team a set of processing tools that can be interconnected at different stages, so that a priori information and expertize are continuously fed to iterative “number crunching” loops. We have analyzed the main problems of SJI, the need for an integration platform and a “proof of concept” of its utilization in two geophysical integration cases. Both, the synthetic example and the application in a real exploration context confirm the effectiveness of our integration approach.

Conflict of Interest

All authors declare no conflicts of interest in this paper.