Structural joint inversion on irregular meshes

SUMMARY

Structural joint inversion of several data sets on an irregular mesh requires appropriate coupling operators. To date, joint inversion algorithms are primarily designed for the use on regular rectilinear grids and impose structural similarity in the direct neighbourhood of a cell only. We introduce a novel scheme for calculating cross-gradient operators based on a correlation model that allows to define the operator size by imposing physical length scales. We demonstrate that the proposed cross-gradient operators are largely decoupled from the discretization of the modelling domain, which is particularly important for irregular meshes where cell sizes vary. Our structural joint inversion algorithm is applied to a synthetic electrical resistivity tomography and ground penetrating radar 3-D cross-well experiment aiming at imaging two anomalous bodies and extracting the parameter distribution of the geostatistical background models. For both tasks, joint inversion produced superior results compared with individual inversions of the two data sets. Finally, we applied structural joint inversion to two field data sets recorded over a karstified limestone area. By including geological a priori information via the correlation-based operators into the joint inversion, we find P-wave velocity and electrical resistivity tomograms that are in accordance with the expected subsurface geology.

Inverse theory, Joint inversion, Tomography

1 INTRODUCTION

Joint inversions have become a popular tool in geophysics to search for models of the subsurface that are in agreement with data acquired with different geophysical methods (Moorkamp et al. 2016). Combining data sets with complementary sensitivity and resolution properties in a joint inversion can mitigate problems related to inversion such as limited resolution and non-uniqueness. Furthermore, models obtained from joint inversion are inherently consistent with multiple data types and can thus improve interpretation and classification as well as help to establish petrophysical relationships between geophysical parameters.

For joint inversions, the property models related to the different data sets need to be coupled. Coupling is either achieved by direct petrophysical relationships or by structural constraints. Petrophysical links prescribe some (often empirical) relation between the physical properties under consideration and have been shown to produce good results, if the petrophysical relationship is reliable (e.g. Heincke et al. 2006; Moorkamp et al. 2011; Lelièvre et al. 2012). However, the petrophysical relations and parameters are often unknown or might only be valid for the specific case under consideration and may fail elsewhere. Compared to petrophysical coupling, structural coupling is a milder and more flexible constraint that only enforces structural similarity between the property models, which makes this approach more widely applicable. The underlying assumption of structural similarity is that the variations in the different geophysical properties originate from the same (geological) structures (Haber & Oldenburg 1997).

Structural joint inversion has been applied to problems over a large range of scales (Moorkamp et al. 2016) and different approaches for structural integration have been proposed in the literature. For example, Günther & Rücker (2006) structurally link two models by changing the regularization for one model according to the structure of the other model. This method has been further refined and applied in several studies (Hellman et al. 2017; Ronczka et al. 2017; Skibbe et al. 2018). Another approach is suggested by Haber & Gazit (2013). They compute a structural measure, called joint total variation, that is based on the norm of the absolute spatial gradient values of the considered models. Yet a different approach for structurally linking two parameter fields is proposed by Lien (2013). In her approach, structure is considered as transition between dominating parameter values and directly inverted for as a common model parameter that represents the structure in the different parameter fields. Thus, additional penalizing terms in the objective function are avoided.

The most popular and widely used structural coupling mechanism is the so-called cross-gradient approach (Gallardo & Meju 2003). Cross-gradient constraints have been applied successfully in joint inversions on a variety of scales and targets, such as in hydrological near-surface studies (e.g. Linde et al. 2006; Doetsch et al. 2010b), targets in mineral exploration (e.g. Lelièvre et al. 2012) or hydrocarbon exploration (e.g. Colombo et al. 2012). The majority of the 2-D applications (e.g. Gallardo 2007; Hu et al. 2009; Colombo et al. 2012) as well as 3-D examples (Linde et al. 2006; Fregoso & Gallardo 2009; Doetsch et al. 2010b; Moorkamp et al. 2011; Shi et al. 2017) use regular grid discretizations. Only Lelièvre et al. (2012) work on irregular tetrahedral (triangular in 2-D) grids. They show an application of structural joint inversion using cross-gradients on irregular meshes in a 2-D synthetic example.

While irregular meshes are superior in representing complicated subsurface structures, topography and complex acquisition geometries (Rücker et al. 2006), the design of gradient and regularization operators is not trivial for irregular discretizations (Lelièvre & Farquharson 2013). A suitable option for regularization on irregular meshes are geostatistical regularization operators (Jordi et al. 2018). Geostatistical regularization is based on a correlation model and allows to incorporate prior knowledge such as length scales over which geological properties are expected to be correlated (correlation lengths). In particular, for geological settings with a preferential direction (e.g. layering), geostatistical regularization can improve the inversion results compared to traditional cell-to-cell smoothness constraints (e.g. Hermans et al. 2012; Jordi et al. 2018).

Concerning the model gradients that are needed to calculate the cross-gradients for structural coupling, Lelièvre & Farquharson (2013) propose to calculate the components of the model gradient at the centre of a particular cell by fitting a linear trend through the model values of that cell and its direct neighbours. For each cell, a small inverse problem is solved, resulting in the sought model gradients. We combine the ideas of Lelièvre & Farquharson (2013) for gradient calculation and those of Jordi et al. (2018) for designing the operator size based on a geostatistical model, in order to develop a flexible and robust way of calculating cross-gradients on unstructured meshes for structural coupling. Instead of only using the nearest neighbours of a cell, we use a correlation-based operator design that includes a larger neighbourhood around a particular cell. The main advantage of the correlation-based operators is their reduced dependency on a particular mesh discretization. Due to the underlying correlation model, they act on a physical length-scale rather than on the length-scale that is prescribed by the cell sizes of a mesh. In particular, for irregular meshes, where cell sizes and distances between neighbouring cells can vary, decoupling of the operators from the mesh is important.

Using these operators, we present a structural joint inversion algorithm that is suitable for 2-D and 3-D joint inversions on irregular meshes. We tested our joint inversion on a 3-D synthetic cross-well example, typical for hydrological near-surface studies. In a second example, joint inversion was applied to 2-D seismic refraction and ERT data that were recorded in a complex karstified area.

2 METHODS

2.1 Joint inversion scheme

The objective function of the joint inverse problem has the following form:

$$\begin{eqnarray*} \Phi = \Phi _d + \lambda \Phi _m + \lambda _{cg}\Phi _{cg} \end{eqnarray*}$$

(1)

where λ and λ_cg are the geostatistical regularization and cross-gradient weights, respectively. Φ_d is the combined data misfit term for all data types

$$\begin{eqnarray*} \Phi _d = \sum _{k=1}^{K}\left\Vert \mathbf {C}_{d,k}^{-0.5}(\mathbf {d}_{k}-\mathbf {F}_k(\mathbf {m}_k) \right\Vert ^2_2, \end{eqnarray*}$$

(2)

where K is the number of data types in the joint inverse problem. We consider here only the case K=2, but extension to more data types is possible (Gallardo 2007; Doetsch et al. 2010b). |$\mathbf {C}_{d,k}$| is the data covariance matrix (assumed to be diagonal), |$\mathbf {d}_{k}$| are the observed data, |$\mathbf {F}_k(\mathbf {m}_k)$| is the calculated forward response of the model |$\mathbf {m}_k$|⁠. The regularization term is defined as

$$\begin{eqnarray*} \Phi _m = \sum _{k=1}^{K}\Vert \mathbf {C}_{m,k}^{-0.5}\Delta \mathbf {m}_k \Vert ^2_2. \end{eqnarray*}$$

(3)

|$\mathbf {C}_{m,k}^{-0.5}$| is the regularization operator for model |$\mathbf {m}_{k}$| and |$\Delta \mathbf {m}_k$| is the difference between the model |$\mathbf {m}_{k}$| and a reference model |$\mathbf {m}^\text{ref}_k$|⁠. The cross-gradients term for two data types (K= 2) and the corresponding models |$\mathbf {m}_1$| and |$\mathbf {m}_2$| is

$$\begin{eqnarray*} \Phi _{cg} = \Vert \mathbf {t}(\Delta \mathbf {m}_1,\Delta \mathbf {m}_2) \Vert ^2_2 , \end{eqnarray*}$$

(4)

where |$\mathbf {t}(\Delta \mathbf {m}_1,\Delta \mathbf {m}_2)$| is the cross-product calculated between |$\Delta \mathbf {m}_1$| and |$\Delta \mathbf {m}_2$|⁠, that is the model updates with respect to the reference model (see below).

For two different data types |$\mathbf {d}_1$| (e.g. apparent resistivity) and |$\mathbf {d}_2$| (e.g. traveltime), the following system of equations is solved in a least-squares sense (adapted from Linde et al. 2006):

$$\begin{eqnarray*} {\begin{bmatrix}\mathbf {C}_d^{-0.5}\mathbf {J}^p \\ \mathbf {C}_m^{-0.5} \\ \lambda _{cg}\mathbf {B}_x^p \\ \lambda _{cg}\mathbf {B}_y^p \\ \lambda _{cg}\mathbf {B}_z^p \end{bmatrix}} {\begin{bmatrix}\Delta \mathbf {m}^{p+1}\\ \end{bmatrix}} = {\begin{bmatrix}\mathbf {C}_d^{-0.5}\skew{3}\hat{\mathbf {d}}^p\\ \mathbf {C}_m^{-0.5} \Delta \mathbf {m}^{p}\\ \lambda _{cg}(\mathbf {B}_x^{p}\Delta \mathbf {m}^{p}-\mathbf {t}_x^{p}) \\ \lambda _{cg}(\mathbf {B}_y^{p}\Delta \mathbf {m}^{p}-\mathbf {t}_y^{p}) \\ \lambda _{cg}(\mathbf {B}_z^{p}\Delta \mathbf {m}^{p}-\mathbf {t}_z^{p}) \end{bmatrix}} , \end{eqnarray*}$$

(5)

where superscript p indicates the iteration number,

$$\begin{eqnarray*} \Delta \mathbf {m}^{p+1} = {\begin{bmatrix}\Delta \mathbf {m}_1^{p+1}\\ \Delta \mathbf {m}_2^{p+1} \end{bmatrix}}, \mathbf {J}^p = {\begin{bmatrix}\mathbf {J}_1^p\\ \mathbf {J}_2^p \end{bmatrix}}, \end{eqnarray*}$$

$$\begin{eqnarray*} \mathbf {\hat{d}}^p = {\begin{bmatrix}\mathbf {d}_1-\mathbf {F}_1(\mathbf {m}^p_1)+\mathbf {J}_1^p\Delta \mathbf {m}_1^{p}\\ \mathbf {d}_2-\mathbf {F}_2(\mathbf {m}^p_2)+\mathbf {J}_2^p\Delta \mathbf {m}_2^{p} \end{bmatrix}}, \end{eqnarray*}$$

$$\begin{eqnarray*} \mathbf {C}_d^{-0.5} = {\begin{bmatrix}w_1\mathbf {C}_{d1}^{-0.5} & \quad 0\\ 0 &\quad w_2\mathbf {C}_{d2}^{-0.5} \end{bmatrix}}, \end{eqnarray*}$$

$$\begin{eqnarray*} \mathbf {C}_m^{-0.5} = {\begin{bmatrix}\lambda ^pw_1\mathbf {C}_{m1}^{-0.5} & \quad 0\\ 0 & \quad \lambda ^pw_2\mathbf {C}_{m2}^{-0.5} \end{bmatrix}}. \end{eqnarray*}$$

|$\Delta \mathbf {m}^{p+1}_1$| (e.g. resistivity) and |$\Delta \mathbf {m}^{p+1}_2$| (e.g. slowness) are the sought model updates that are used to iteratively update the model vectors with respect to the reference models |$\mathbf {m}_{1,2}^{\rm {ref}}$| (Doetsch et al. 2010a):

$$\begin{eqnarray*} {\begin{bmatrix}\mathbf {m}_1^{p+1}\\ \mathbf {m}_2^{p+1} \end{bmatrix}} = {\begin{bmatrix}\mathbf {m}_1^{{\rm {ref}}}\\ \mathbf {m}_2^{{\rm {ref}}} \end{bmatrix}}+ {\begin{bmatrix}\Delta \mathbf {m}_1^{p+1}\\ \Delta \mathbf {m}_2^{p+1} \end{bmatrix}}. \end{eqnarray*}$$

(6)

|$\mathbf {J}^p_1$| and |$\mathbf {J}^p_2$| are the Jacobian matrices containing the sensitivities of a change in data |$\mathbf {d}_{1,2}$| related to a change in model |$\mathbf {m}_{1,2}$|⁠, w₁ and w₂ are weighting factors (see below), λ^p is the overall regularization weight that is allowed to change for each iteration and λ_cg is the cross-gradient weight. |$\mathbf {B}^{p}_{x,y,z}$| are the Jacobians of the cross-gradient functions |$\mathbf {t}^{p}_{x,y,z}$| with respect to changes in the respective models and respective directions. These sensitivities, for example, |$\mathbf {B}^{p}_{x}$|⁠, are used to estimate the cross-gradient function |$\mathbf {t}_x^{p+1}$| at the present iteration by a first order Taylor expansion around the cross-gradient function |$\mathbf {t}_x^{p}$| of the previous iteration (Linde et al. 2008):

$$\begin{eqnarray*} \mathbf {t}_x^{p+1} \cong \mathbf {t}_x^{p} + \mathbf {B}_x^{p}\Delta \mathbf {m}^{p}. \end{eqnarray*}$$

(7)

Weighting the different data sets is a key issue in joint inversion (Doetsch et al. 2010b; Moorkamp et al. 2016; Steklova & Haber 2017). We propose a weighting scheme that uses the data fit measure χ² = Φ_d/n_d, where n_d is the number of data. The weights w₁ and w₂ are calculated as follows:

$$\begin{eqnarray*} w_1 = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1.0 & \quad \text{if } \chi ^2_1 \gt \chi ^2_2\\ \chi ^2_1/\chi ^2_2 & \quad \text{if } \chi ^2_1 \lt \chi ^2_2 \end{array}\right. \end{eqnarray*}$$

and

$$\begin{eqnarray*} w_2 = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1.0 & \quad \text{if } \chi ^2_1 \lt \chi ^2_2\\ \chi ^2_2/\chi ^2_1 & \quad \text{if } \chi ^2_1 \gt \chi ^2_2 \end{array}\right. \end{eqnarray*}$$

(8)

This weighting scheme should balance the convergence of the individual data misfits in eq. (5).

2.2 Regularization operators

We use the geostatistical regularization operators described by Jordi et al. (2018). These operators can be thought of as superposition of damping and smoothing constraints (Maurer et al. 1998). Unlike standard regularization operators that are defined using the nearest neighbours of a cell only, geostatistical operators include a larger neighbourhood around each cell, where a covariance matrix |$\mathbf {C}_m$| contains the correlation information that defines the extent of this neighbourhood. We use the exponential correlation model for which the covariance matrix is found as:

$$\begin{eqnarray*} \mathbf {C}_m=\sigma ^{2}\exp {\left(-3\sqrt{\left(\frac{\mathbf {H}_x}{I_x}\right)^2+\left(\frac{\mathbf {H}_y}{I_y}\right)^2+\left(\frac{\mathbf {H}_z}{I_z}\right)^2}\right)} , \end{eqnarray*}$$

(9)

where the matrices |$\mathbf {H}_x$|⁠, |$\mathbf {H}_y$| and |$\mathbf {H}_z$| contain the lag distances, that is the distances of every cell centre to all other cell centres. σ² is the variance of the model parameters and I_x, I_y and I_z are the correlation lengths in the x, y and z directions, respectively (e.g. Kitanidis 1997). The correlation length defines the distance over which cells are considered to be correlated, that is, the distance for which the correlation function decays to 5 per cent. The regularization operators |$\mathbf {C}^{-0.5}_m$| used in eq. (5) are found via an eigendecomposition of the covariance matrix |$\mathbf {C}_m$| (Jordi et al. 2018) and are imposed instead of conventional regularization (damping or smoothing).

2.3 Gradients, cross-gradients and cross-gradient sensitivities

In order to calculate the cross-gradients in 3-D, the x-, y- and z-components of the local model gradient vector are required. Cross-gradient calculations designed for regular grids (e.g. Gallardo & Meju 2003; Linde et al. 2008; Doetsch et al. 2010a) calculate differences between adjacent cells to determine the gradient components. We do not rely on the direct calculation of model differences and rather adapt the approach of Lelièvre & Farquharson (2013) for the calculation of the local spatial gradients in the model updates |$\Delta \mathbf {m}$|⁠. They calculate the model gradient at the location of a cell by fitting a linear trend through the model parameters of that cell and its neighbouring cells. This method avoids patchwork problems that could appear otherwise (Lelièvre & Farquharson 2013). We generalize this approach and use all cells that lie within a predefined distance (neighbourhood) for the gradient calculation. The equation of the linear trend for a cell i with model update value Δm_i in a 3-D model is

$$\begin{eqnarray*} \Delta m_{i}(x_i,y_i,z_i) = ax_i + by_i + cz_i + d , \end{eqnarray*}$$

(10)

where the linear trend parameters a, b and c represent the x-, y- and z-component of the spatial gradient |$\mathbf {g}_i = (a,b,c)^T$|⁠, respectively. For a neighbourhood of n model cells, a system of equations of the form |$\mathbf {X}_i \mathbf {q}_i = \Delta \mathbf {\hat{m}}_{i}$| can be formulated

$$\begin{eqnarray*} {\begin{bmatrix}x_i & \quad y_i & \quad z_i & \quad 1 \\ x_1 & \quad y_1 & \quad z_1 & \quad 1 \\ x_2 & \quad y_2 & \quad z_2 & \quad 1 \\ \vdots & \quad \vdots & \quad \vdots & \quad \vdots \\ x_n & \quad y_n & \quad z_n & \quad 1 \end{bmatrix}} {\begin{bmatrix}a \\ b \\ c \\ d \end{bmatrix}} = {\begin{bmatrix}\Delta m_i \\ \Delta m_1 \\ \Delta m_2 \\ \vdots \\ \Delta m_n \end{bmatrix}} , \end{eqnarray*}$$

(11)

where the vector |$\mathbf {q}_i$| contains the gradient components and matrix |$\mathbf {X}_i \in \mathbb {R}^{n \times 4}$| contains the coordinates of the n model cells that contribute to the calculation of the gradient. |$\Delta \mathbf {\hat{m}}_{i}$| is a vector holding the respective model update values of these cells. This system is overdetermined if at least the cell for which the gradient is calculated (cell i in eq. (11)) and its direct neighbours are used to fit the linear trend. In case of a 3-D tetrahedral mesh, cell i and the four neighbouring cells [n = 4 in eq. (11)] are needed to calculate the spatial gradient. The gradient components can then be found through a least squares inversion (Lelièvre & Farquharson 2013):

$$\begin{eqnarray*} \mathbf {q}_i = \mathbf {X}_i^{-1}\Delta \mathbf {\hat{m}}_{i} . \end{eqnarray*}$$

(12)

We introduce a weighting operator |$\mathbf {W}_i \in \mathbb {R}^{n \times n}$| that is based on a correlation model. As for the regularization operator, we use an exponential correlation model (eq. 9), which is decreasing the weight for the gradient calculation of cells with increasing distance from cell i. The cells contributing to the gradient calculation (neighbourhood) are restricted to lie within an ellipsoid (ellipse in 2-D) with semi-axes equal to one correlation length in the respective direction. The operator |$\mathbf {W}_i$| is a diagonal matrix containing the truncated correlation model of cell i (row i of |$\mathbf {C}_m$|⁠). Note that both |$\mathbf {X}_i$| and |$\mathbf {W}_i$| are purely dependent on the mesh geometry and thus, if the mesh is not changed during the inversion process, they need to be determined only once at the beginning of the inversion. For each cell an operator of the following form is determined:

$$\begin{eqnarray*} \mathbf {\hat{X}}_i = \left( \mathbf {X}_i^T\mathbf {W}_i^T\mathbf {X}_i \right)^{-1}\mathbf {X}_i^T\mathbf {W}_i^T\mathbf {W}_i , \end{eqnarray*}$$

(13)

where |$\mathbf {\hat{X}}_i \in \mathbb {R}^{3 \times n}$| is reduced to the columns needed for the calculation of the three gradient components. The gradient vector |$\mathbf {g}(\Delta \mathbf {\hat{m}}_{i})$| at the location of a cell i is then found by

$$\begin{eqnarray*} \mathbf {g}(\Delta \mathbf {\hat{m}}_{i}) = \mathbf {\hat{X}}_i\Delta \mathbf {\hat{m}}_{i}. \end{eqnarray*}$$

(14)

The cross-gradient for cell i and two model updates |$\Delta \mathbf {m}_{1}$| and |$\Delta \mathbf {m}_{2}$| is then:

$$\begin{eqnarray*} \mathbf {t}_i(x,y,z) = \mathbf {g}(\Delta \mathbf {\hat{m}}_{1,i}) \times \mathbf {g}(\Delta \mathbf {\hat{m}}_{2,i}). \end{eqnarray*}$$

(15)

Also for the calculation of the cross-gradient sensitivity matrices |$\mathbf {B}_{x,y,z}$| in eq. (5), we make use of the geometry-dependent operator |$\mathbf {\hat{X}}_i$|⁠. The cross-gradient sensitivity for cell i with respect to a change in model value m_{1, j} are calculated as follows:

$$\begin{eqnarray*} \mathbf {B}_{1,i} &=& \dfrac{\partial t_i}{\partial m_{1,j}} = \dfrac{\partial (\mathbf {g}(\Delta \mathbf {\hat{m}}_{1,i}) \times \mathbf {g}(\Delta \mathbf {\hat{m}}_{2,i}))}{\partial m_{1,j}} \nonumber \\ &=& \mathbf {g}(\Delta \mathbf {\hat{m}}_{1,i}) \times \dfrac{\partial (\mathbf {g}(\Delta \mathbf {\hat{m}}_{2,i}))}{\partial m_{1,j}} + \dfrac{\partial (\mathbf {g}(\Delta \mathbf {\hat{m}}_{1,i}))}{\partial m_{1,j}} \times \mathbf {g}(\Delta \mathbf {\hat{m}}_{2,i})\nonumber \\ &=& \mathbf {\hat{X}}_{ij} \times \mathbf {g}(\Delta \mathbf {\hat{m}}_{2,i}) , \end{eqnarray*}$$

(16)

where |$\mathbf {\hat{X}}_{ij}$| represents the jth row of the operator matrix |$\mathbf {\hat{X}}_{i}$|⁠. Similarly, the sensitivity to a change in model value m_{2, j} is calculated:

$$\begin{eqnarray*} \dfrac{\partial t_i}{\partial m_{2,j}} = \mathbf {g}(\Delta \mathbf {\hat{m}}_{1,i}) \times \mathbf {\hat{X}}_{ij}. \end{eqnarray*}$$

(17)

3 EXAMPLE OF CORRELATION-BASED CROSS-GRADIENT CALCULATION

We examine the calculation of the cross-gradients using the correlation-based approach described above on two synthetic 2-D models (Fig. 1a). The models are composed of two sinusoidal variations. A horizontal and a vertical sinusoidal of wavelengths 80 and 20 m, respectively. The models mimic a setting where a preferential direction, for example, layering, is present.

Figure 1.

(a) Synthetic model m₁. (a) Synthetic model m₂. (c) Box plot of differences between analytic and cross-gradients calculated from footprint operator for discretization of m₁ and m₂ with maximum allowed cell sizes of 2, 4 and 8 m². (d) Differences between analytic and cross-gradients calculated from nearest neighbours operator. The upper and lower limits of the black box are the first and third quartile, defining the interquartile-range (IQR). The orange line within the box represents the second quartile (median). The black bars indicate the lowest datum and highest datum that are 1.5 IQR within the lower and the upper quartile, respectively. Outliers are marked as black circles.

For both models, analytic gradient fields can be determined and used to calculate analytic cross-gradients. The cross-gradients are evaluated analytically at the cell centres and compared to cross-gradient fields calculated either with a predefined operator footprint or only based on the nearest neighbours of a cell. Moreover, cross-gradients are calculated for different mesh discretizations, where the maximum cell sizes are restricted to 2, 4 and 8 m². In Figs 1(a) and (b), only the discretization with maximum cell sizes of 4 m² is shown. The box plot in Fig. 1(c) shows the discrepancy between the analytic cross-gradient field and the cross-gradients calculated with an operator with correlation lengths of I_x = 8 m and I_z = 2 m. Similarly, for an operator comprising only the nearest neighbours, the differences to the analytic cross-gradients are shown in Fig. 1(d). Below the two box plots, the typical operator sizes are shown for the meshes with different maximum cell sizes.

First, we observe in the plots of the operators that the correlation-based operators vary less in size compared to the operators with only the nearest neighbours. The latter scale with increasing cell sizes, while the size of the former remains approximately constant for the different cell sizes. This is due to the fact that the correlation-based operators are designed to have a mesh-independent size. From Figs 1(c) and (d) we observe that compared to the nearest neighbour case, the calculated cross-gradients are closer to the analytic cross-gradients for the correlation-based operators. Furthermore, the deviation of the calculated cross-gradients from the analytic ones shows a stronger increases with increasing cell size for the nearest neighbour case. The interquartile range (IQR, see Table 1) of the calculated deviation increases for the footprint operator by 14.5 per cent when the cell size is increased from 2  to 4 m². For the nearest neighbour operator the IQR increases by 36.5 per cent. An increase in cell size from 2  to 8 m² results in a 67.6 and 91.1 per cent larger IQR for footprint and nearest neighbour operators, respectively. Also when comparing the values of ϕ_cg we observe that the relative difference to the value of the analytic cross-gradients only slightly varies for the footprint operators but exhibits large changes for the nearest neighbour operators (see Table 1). These relative differences are 5.9 per cent (2 m²), 5.5 per cent (4 m²) and 5.0 per cent (8 m²) in case of the footprint operators and – 17.9 per cent (2 m²), – 34.7 per cent (4 m²) and – 106.8 per cent (8 m²) for the nearest neighbour operators. We interpret this behaviour as an indication that the correlation-based approach is less dependent on the mesh discretization.

4 JOINT INVERSION OF SYNTHETIC 3-D CROSS-HOLE DATA

4.1 Synthetic model

Our 3-D joint inversion algorithm was tested with a synthetic cross-hole electrical resistivity tomography (ERT) and ground-penetrating radar (GPR) traveltime tomography scenario that is inspired by a field study presented by Doetsch et al. (2010b). The two models consist of layered background variations and two homogeneous anomalies (Fig. 2). The layered background has been generated using geostatistical simulation of porosity ϕ values that are typically found in sandy aquifer environments. The porosity field was then used to calculate the resistivity model ρ via Archie’s law (Archie 1942), assuming full saturation:

$$\begin{eqnarray*} \rho = \phi ^{-m}\rho _w \end{eqnarray*}$$

(18)

with the cementation exponent m =  1.8 and the resistivity of water ρ_w  =  60 Ωm. The GPR-velocity values were calculated using the complex refraction index model (CRIM, Lesmes & Friedman 2005):

$$\begin{eqnarray*} \sqrt{k_f}=\phi \sqrt{k_w} + (1-\phi )\sqrt{k_m} \end{eqnarray*}$$

(19)

where k_w  =  81 and k_m = 3 are the relative permittivities of water and the matrix material. They represent typical values for a sandy aquifer. The total relative permittivity k_f is then used to calculate the GPR-velocities v_r:

$$\begin{eqnarray*} v_r=c/\sqrt{k_f}, \end{eqnarray*}$$

(20)

where c is the speed of light.

Figure 2.

(a) Synthetic resistivity model. (b) Synthetic GPR velocity model. The black spheres are electrode and GPR antenna positions, grey spheres are GPR antenna positions only.

Two anomalies were superimposed on the layered background models. A sphere centred at (2.5, 2.5 and −1.85 m) with radius 1.2 m, and an ellipsoid with midpoint (2.5, 2.75 and −5.1 m) and semi-major axes lengths of 2.0, 1.5 and 1.0 m in x-, y- and z-direction, respectively. The sphere represents a high velocity (85 m μs^–1) and high resistivity anomaly (370 Ωm), whereas the ellipsoid is a low velocity (55 m μs^–1) and low resistivity (120 Ωm) anomaly.

The acquisition setup and measurement schemes for the synthetic GPR and ERT data were set up as follows. Four boreholes are placed at the corners of a 5 m × 5 m square (Fig. 2). Each borehole contains nine electrodes evenly spaced with 0.7 m vertical separation, starting from z = –0.7 m down to z = −6.3 m (black spheres in Fig. 2). For the GPR measurements, a total of 20 GPR transmitter and receiver positions spaced at 0.35 m are simulated in each borehole (black and grey spheres in Fig. 2). The topmost GPR antenna position in each borehole is at z = 0.0 m and the deepest at z = –6.65 m. For ERT we considered both AB-MN and AM-BN recording schemes (see Doetsch et al. 2010b, for details) with electrodes located in two different boreholes. All possible pairs of boreholes have been used to model 5000 ERT recordings that were contaminated with 3 per cent Gaussian noise. For the GPR first arrival traveltimes, all possible cross-hole combinations have been simulated, resulting in 4800 traveltimes. The traveltime data were contaminated with absolute Gaussian noise (standard deviation of 1 ns).

4.2 Results and model assessment

We compare the joint inversion results to (i) those obtained by individually inverting each data set and (ii) to the true models (Fig. 2). For the individual inversions, we used the same framework as described above for joint inversion (eq. 5) but using only one method, for example those parts related to ERT and omitting all terms related to the cross-gradient coupling. For ERT, the model responses and sensitivities were calculated using the finite-element algorithms described in Rücker et al. (2006) and Günther & Rücker (2006), which are implemented in the software package pyGIMLi (Rücker et al. 2017). Similarly, for GPR-traveltime tomography, the inversion scheme was only applied to the traveltime data. The model responses were calculated using a Shortest Path algorithm (Dijkstra 1959), that is also embedded in pyGIMLi. In order to obtain an estimation of the correlation parameters for geostatistical regularization, a variogram analysis in vertical direction of the true background resistivity and GPR-slowness models was performed (anomalies not included). This analysis mimics a typical case where limited prior logging information from the existing boreholes is used to infer subsurface properties. The variogram analysis yielded an estimated variance for the logarithms of the resistivity (Ωm) and GPR-slowness (μs m^–1), which are the properties we invert for, of 0.03 and 0.006, respectively. The vertical correlation length was determined to be 2 m. For the horizontal correlation lengths we used 5 m which is our best estimation if we correlate the layers across the different boreholes, not knowing about the presence of the anomalies in the centre of the domain that intersect these layers.

The results from individually inverting the synthetic ERT and traveltime data sets are shown in Figs 3(a) and (b). In Fig. 3(a) we observe that the layered background resistivity model is well recovered. However, only the low resistivity anomaly in the lower part of the domain is visible. Due to decreasing sensitivity away from the boreholes (Day-Lewis et al. 2005) and generally low sensitivity for high resistivity bodies, the high resistivity anomaly is almost absent in the ERT tomogram. In the low sensitivity areas the geostatistical regularization enforces the layering according to the correlation lengths I_x = I_y= 5 m and I_z = 2 m that were used (Jordi et al. 2018). Within 6 iterations and a regularization weight of λ = 3 the inversion converged to a data fit of χ² = 1.0, explaining the data to the error level. The regularization operator for GPR traveltime tomography (Fig. 3b) has the same correlation lengths as for the ERT inversion but the estimated variance is 0.006. Setting λ = 4, four iterations were needed to reach a χ² of 1.0. Using traveltime tomography, both anomalies as well as the layered background are imaged. Around the high velocity anomaly, smearing of high velocity values due to the geostatistical regularization is observed.

Figure 3.

(a) Resistivity tomogram from individual inversion. (b) GPR velocity tomogram from individual inversion. (c) Resistivity tomogram from joint inversion. d. GPR velocity tomogram from joint inversion. The dashed white lines indicates the true positions of the two anomalies.

The results of jointly inverting the two data sets are shown in Figs 3(c) (ERT) and (d) (GPR). The weights for the individual data sets w₁ (ERT) and w₂ (GPR) were determined via the data fit χ²-criterion (eq. 8) and to enforce a high degree of structural similarity, the cross-gradient weight was set to λ_cg = 8×10⁴ after testing values in the range of 1×10⁴ to 1×10⁵. The differences between inversion results for different λ_cg are marginal. With the chosen λ_cg good convergence and minimization of the cross-gradients was achieved. The regularization operator used for joint inversion was the same as for the individual inversions. The overall regularization weight λ was set to 120 and automatically decreased by 20 per cent after each iteration down to λ of 4.2 for iteration 16. We found that decreasing a high initial λ between 100–200 yields good inversion results. The inversion is terminated if either of the data is fitted to a χ² of 1 or a pre-defined maximum number of iterations is reached. Within 16 iterations, both, ERT and GPR traveltime data were fitted to a χ² of 1.0. The two convergence curves of the data fit χ² for 20 iterations are shown in Fig. 4. Similar convergence behaviour for both data sets is observed. Along with the data fit, the cross-gradient term of the objective function (eq. 4) is displayed as |$\Phi _{cg}*\lambda _{cg}^2$|⁠. Starting from the two homogeneous initial models where |$\Phi _{cg}*\lambda _{cg}^2$| is zero (iteration 0), structure is introduced into the models which is reflected by an increase of |$\Phi _{cg}*\lambda _{cg}^2$|⁠, up to a maximum structural dissimilarity at iteration 8. From iteration 9 onwards, |$\Phi _{cg}*\lambda _{cg}^2$| is progressively decreased. At this point, the previously fast decrease in data misfit χ² is slowed down and structural similarity between the models is established through the minimization of the cross-gradients. For our preferred models at iteration 16 (red circles in Fig. 4) |$\Phi _{cg}*\lambda _{cg}^2$| is 84. Compared to the previous iterations, only a small decrease in this value is observed for the following iterations (55 at iteration 20).

Figure 4.

Convergence curves of the data misfit χ² and the cross-gradient term of the objective function Φ_cg (eq. 4) for the joint inversion of the 3-D ERT and traveltime data. The red circles indicate the respective values of χ² and |$\Phi _{cg}*\lambda _{cg}^2$| for the preferred models after 16 iterations (Figs 3c and d).

Comparing the joint inversion results in Figs 3(c) and (d) to the results of the individual inversion, we observe that there are only small differences between the two GPR-velocity tomograms (Figs 3b and d). The high velocity anomaly is slightly better localized, when the joint inversion is used. The differences between the resistivity tomograms are more pronounced. Using joint inversion, the high resistivity anomaly is better visible in the ERT tomogram and the layers towards the centre of the domain are better imaged (Fig. 3c).

These observations are underlined by the differences calculated between the true model and the respective inversion results (Fig. 5). For the calculation of the model differences, the property values of the true models were interpolated on the inversion mesh. While for the difference of the individual ERT inversion in Fig. 5(a) the misfit for the high resistivity anomaly is 150 Ωm, it is decreased to values below 75 Ωm for joint inversion (Fig. 5c). However, enforcing structural similarity in the joint inversion introduced the artefacts of high velocity smearing observed for traveltime tomography into the ERT tomogram. Around the high resistivity anomaly, a halo of negative difference values is visible in Fig. 5(c), which is not present for the individual ERT inversion (Fig. 5a). As already observed, there are no significant differences between individual and joint inversion results for the GPR-velocity tomograms (Figs 5b and d).

Figure 5.

Difference between true and inverted models. (a) Resistivity tomogram from individual inversion. (b) GPR-velocity tomogram from individual inversion. (c) Resistivity tomogram from joint inversion. d GPR-velocity tomogram from joint inversion.

Although the differences in the tomograms are small, the joint inversion has a significant impact on the petrophysical distribution of the two parameter fields. This is visualized in the scatter plot of log-resistivity versus log-GPR-slowness (Fig. 6a). The black plus signs represent the distribution of the true models, orange and blue points are the distributions for the results of joint and individual inversions, respectively. Neither joint, nor individual inversions are able to capture the true resistivity and slowness values of the anomalies, which are the isolated black points in the top left corner (high resistivity and high velocity) and bottom right corner (low resistivity and low velocity) of the scatter plot. The distribution obtained from the joint inversion models (orange) follows the linear trend of the true background model (black) for low resistivities and velocities but starts to deviate from the true trend for high resistivities and velocities towards the values of the high resistivity and high velocity anomaly. The distribution from the individual inversion models scatters around the true distribution. The large scatter and the missing relationship between the two properties makes it difficult to determine a petrophysical relationship from the individual inversion result. The joint inversion result is clearly superior in this respect, as a direct relationship between resistivity and velocity is evident in the scatter plot.

Figure 6.

(a) Scatter plot of the logarithm of resistivity ρ versus the logarithm of slowness 1/v_r. Normalized histogram of the porosities calculated from (b) individual inversion and (c) joint inversion results of GPR (blue) and ERT (orange). The black contour shows the true porosity distribution. The count is normalized by the number of observations times the bin width.

Besides a relationship between the geophysical parameter fields, a translation of the geophysical properties to other properties such as porosity is often of interest. Using Archie’s law (eq. 18) and CRIM (eqs 19 and 20) with the parameters used to calculate the true geophysical property fields (Fig. 2) from the modelled porosity field, we back-calculated the porosity from the individual and joint inversion results. In Fig. 6(b) the normalized porosity distributions calculated from the resistivity (orange) and GPR velocity (blue) results of individual inversions are displayed. In Fig. 6(c) the distributions based on the joint inversion results are shown. The black line represents the true porosity distribution. Note, the true distribution also shows the porosity of the anomalies calculated from GPR velocities at 0.25 and 0.56, which were not modelled by the initial porosity field. As already observed on the scatter plot, neither of the models reaches the maximum or minimum values of the anomalies. However, the two porosity distributions obtained from the joint inversion results Fig. 6(c) are very similar to each other and match the true distribution better compared to the distributions calculated from the models of the individual inversions Fig. 6(b).

5 JOINT INVERSION OF 2-D SURFACE FIELD DATA

We applied our joint inversion approach to two resistivity and traveltime field data sets that were recorded along the same profile in the Swiss Jura mountains for studying karstified Malm limestone units (Schmelzbach et al. 2015). The ERT data set consists of 16 700 measurements recorded with dipole-dipole and multigradient measurement schemes. The estimated error from reciprocal measurements is 3 per cent. The nominal spacing of the 296 electrodes is 2.5 m, covering a profile of 720 m. On the same profile, a seismic survey was carried out. A total of 127 source points (sledgehammer) spaced at 5 m and 296 geophones nominally spaced at 2.5 m were used. The survey resulted in a total of 14 200 first break picks up to a maximum source–receiver offset of 260 m. Based on the picking accuracy of the first P-wave arrivals and a reciprocity analysis (Schmelzbach et al. 2008) a rough estimation of the uncertainty in the traveltime data of 5 per cent was obtained. For the regularization operators we use the parameters that were found using prior information from geological maps and digital elevation models as described in Jordi et al. (2018). The vertical and horizontal correlation lengths were set to 15 and 50 m, respectively. Additionally, the operator was rotated by 12° from the horizontal, which corresponds to the predominant dipping angle of the geological layers. The variance was set to 1/9.

The results from individually inverting the ERT and traveltime data are shown in Figs 7(a) and b. For the ERT inversion, a data fit of χ² = 1.0 was obtained after five iterations. For the velocity model in Fig. 7(b) the data fit is χ² = 1.2. Although it is possible to fit the data to a χ²-value closer to 1, the inversion was terminated after four iterations in order to have a comparable data fit to the joint inversion (see below). We observed only small differences between models with a data fit closer to 1 and the model shown in Fig. 7(b).

Figure 7.

(a) Result from individually inverting the ERT data. Areas of lower sensitivity are progressively faded out with decreasing sensitivity. (b) Result from inverting the seismic traveltime data. Only cells that are covered by rays are displayed. For a better comparison to the result from joint inversion, the tomogram is displayed with the ray coverage of the joint inversion and deeper cells that are covered by rays from individual inversion are semi transparent. ERT tomogram (c) and P-wave velocity tomogram (d) from joint inversion. The black arrow marks the position of a know subvertical fault structure.

The tomograms in Figs 7(a) and (b) share some common features, as already observed in Jordi et al. (2018). The most prominent structure is a high resistivity (ρ > 2000 Ωm) and high P-wave velocity (v_p > 3000 m s^–1) zone between x= 80 and 350 m. Both models show dipping structures for x> 350 m. At around x = 510 m, a zone of low resistivity (ρ <1000 Ωm) and low velocity (v_p <2000 m s^–1) intersects the predominantly layered features in this area.

The structures observed in the tomograms from individual inversions indicate that structural similarity between the two parameter fields exists, and thus, we applied our structural joint inversion algorithm to the two data sets. The results of this joint inversion are shown in Figs 7(c) and (d). The same regularization operators as in the individual inversions were applied to the joint inversion. The overall trade-off parameter between model roughness and data fit was set to λ=140 and decreased by 25 per cent after each iteration but λ was not allowed to fall below 10. The cross-gradient weight λ_cg was set to 4300. The data fit obtained for the ERT and traveltime data is χ²=1.0 and χ²=1.3, respectively. The weighting scheme according to eq. (8) was imposed to balance the contribution of the two methods in the joint inversion. The convergence of the data fit values as well as the cross-gradient part of the objective function are shown in Fig. 8. The inversion was stopped at iteration 45 where the ERT data was fitted to the error level (χ²=1.0). No significant decrease in Φ_cg is observed for iteration numbers larger than 45.

Figure 8.

Convergence curves of the data misfit χ² and the cross-gradient term of the objective function Φ_cg (eq. 4) for the joint inversion of the 2-D ERT and traveltime field data. The red circles indicate the respective values of χ² and |$\Phi _{cg}*\lambda _{cg}^2$| for the preferred models after 45 iterations (Figs 7c and d).

Main features, such as the high resistivity and high velocity body between x = 80 and 350 m, observed in the results from individual inversions are preserved when joint inversion is used (cf. Figs 7c and d). The layers for x> 350 m are enhanced on the tomograms obtained from the joint inversion. Also, the low velocity and low resistivity zone at x = 510 m is clearly visible when joint inversion is applied. The main difference between individually and jointly inverting the data sets are the pronounced low velocity zones that are imaged when using joint inversion (Fig. 7d). On the southeast end of the profile, a low velocity structure (v_p <1000  m s^–1), dipping towards northwest, is observed. At the same location in the P-wave velocity tomogram from individual inversion (Fig. 7b), higher velocity values are present. Also in the area where dipping layers are observed (x > 350 m), joint inversion images structures with lower velocities compared to the individual inversion. For the ERT result, we observe more continuity of the layered structures for x > 350 m when joint inversion is used (Fig. 7c). In particular, a thin layer of higher resistivity values (ρ > 2500 Ωm) overlain by a layer of lower resistivities (ρ <2000 Ωm) between x = 350  and 450 m appears as a continuous structure in the ERT tomogram in Fig. 7(c).

5.1 Interpretation

Our interpretation is based on the results from joint inversion in Figs 7(c) and (d) which we prefer over the results obtained from individual inversions because we consider these results to be in better accordance with the knowledge about regional geology and tectonic features. The predominant layered and karstified limestones that are very typical for the Jura mountains are evidently imaged using joint inversion. In particular, layers of low velocity attributed to karstification are imaged using joint inversion. Nevertheless, main geological features such as a known subvertical fault zone indicated by low resistivity and low velocity values at x = 510 m (black arrow in Fig. 7) or the pronounced high resistivity/velocity body between x = 80  and 350 m are observed on all tomograms. To guide our interpretation, we manually classified velocity and resistivity values in four clusters, displayed in Fig. 9. In the scatter plots from individual inversion (Fig. 9a) and joint inversion (Fig. 9c) the clusters are coloured as follows: red represents resistivity values ρ >2000 Ωm, yellow ρ ≤2000 Ωm and velocities v_p >2800  m s^–1, light green 500>ρ ≤2000 Ωm and v_p ≤2800  m s^–1 and blue ρ ≤500 Ωm and v_p ≤2800  m s^–1. In Figs 9(b) and (d) the subsurface models of individual and joint inversions are shown in the same colour scheme, respectively. The blue colours are interpreted to be the soil cover and heavily weathered areas where conductive clay material might be present. Red areas are accounted to mostly intact limestone. The light green and yellow patches are interpreted as karstified limestone. Where the higher velocities of the yellow areas indicate less karstification compared to areas with lower velocities (light green). Based on the joint inversion results, we interpret in Fig. 9(d) karstified structures reaching to greater depths compared to the individual inversions (Fig. 9b). At the south-east end of the profile, between x = 280 and 400 m and towards the northwest end of the profile pronounced lightgreen and blue regions are observed. Also the fault zone at x= 510 m appears in blue and lightgreen colours in Fig. 9(d).

Figure 9.

Scatter plots of velocity values versus resistivity values of the models obtained from individual (a) and joint inversions (c). Note only cells that are covered by rays are considered. Four regions were manually set. Red: ρ >2000 Ωm; yellow: ρ ≤2000 Ωm and v_p >2800 m s^–1; light green: 500>ρ ≤2000 Ωm and v_p ≤2800 m s^–1; blue: ρ ≤500 Ωm and v_p ≤2800 m s^–1. Panels (b) and (d) subsurface models colour coded according to scatter plot for individual and joint inversions, respectively.

6 DISCUSSION

To date, most structural joint inversion algorithms rely on a regular, rectilinear discretization of the model domain. We extended a scheme of Lelièvre & Farquharson (2013) that was designed for the calculation of gradient (regularization) operators for the calculation of cross-gradient operators on irregular triangular and tetrahedral meshes. Unlike any other cross-gradient joint inversion, our algorithm approach offers a way to impose structural similarity between different property fields on a defined length scale. In particular for irregular meshes, where cell sizes vary and cells are oriented randomly, our approach can stabilize the cross-gradient calculation and reduce the dependency of the cross-gradient operators on a particular mesh discretization. If cell sizes are large in relation to the correlation lengths, the operator might include only a few additional cells compared to a nearest neighbour operator, downgrading the performance of the footprint cross-gradient operator. However, in any of the studied cases where we compared the calculated to analytic cross-gradients (Fig. 1), the footprint operator performed better than the nearest neighbour operator.

Although we have not fully optimized the code, we provide here a brief comparison on computational cost and memory requirements of our joint inversion approach. All computations were run on an 8-core 3.5 GHz machine with 64 GB RAM. For the same 3-D cross-hole setting as described above, a refined mesh comprising 20459 cells was generated. The memory requirements and average computation times (over 20 iterations) for inversions with three differently sized isotropic cross-gradient operators are summarized in Table 2. With increasing operator size the number of cells per operator is growing and the LSQR as well as the calculation of the cross-gradient sensitivities become more time consuming. For the smaller operators (correlation lengths of 0.5 and 1 m) most of the iteration time is spent for those parts of the inversion loop that are not dependent on operator size, like calculating the Jacobian matrices (on average 3.2 min) or the forward responses (on average 1.2 min). While time consumption, particularly for the parallelized cross-gradient sensitivity calculations, can be reduced by running the inversion on a cluster with more nodes, the high demand in memory of the current implementation limits the capability of the 3-D joint inversion to near-surface applications with a few ten thousand model parameters. For the small operators the maximum memory demand (28.8 GB) is dictated by the geostatistical regularization operator calculated before the inversion (13.2 min). For the large operator (correlation lengths of 2 m) the number of cells comprising the operator is increased and thus the memory requirement for the cross-gradient sensitivities is significantly higher (41.8 GB). Future work will be dedicated to optimize the computation, particularly the cross-gradients sensitivity calculation.

Table 1.

Interquartile ranges (IQR) and values of Φ_cg for different cell sizes and operators.

Mesh	IQRFP	IQRNN	Φcg, A	Φcg, FP	Φcg, NN
2m2	24.9	239.9	8.336e+8	4.288e+8	1.998e+8
4m2	28.5	327.5	7.848e+8	4.042e+8	1.897e+8
8m2	41.7	458.4	9.829e+8	5.777e+8	4.131e+8

Mesh indicates the max. cell size allowed in mesh. IQR for footprint (FP) and nearest neighbour (NN) operators. Φ_cg of analytic (A), footprint and nearest neighbour cross-gradients.

Table 2.

Average time and memory requirements per iteration for 3-D joint inversion with differently sized operators.

Operator size	avg. cells	avg. time	time LSQR	time CG-sens.	max. memory
0.5/0.5/0.5	47	10.5	1.2	0.4	28.8
1.0/1.0/1.0	308	13.9	2.3	2.4	28.8
2.0/2.0/2.0	1914	37.6	11.3	15.2	41.8

Operator size: Correlation lengths I_x/I_y/I_z in [m]. avg. cells: Average number of cells per operator. avg. time: Average time in [min] per iteration. time LSQR: Average time in [min] for LSQR. time CG-sens.: Average time in [min] for cross-gradients sensitivity calculation. max. memory: Maximum memory in [GB] required.

In this study we defined the cross-gradient operators on the same scale as the regularization, which are scaled according to prior knowledge. There might exist scenarios where the scale of the regularization and structural similarity are not coinciding. Future work will be required to understand, how structural similarity on different length scales and the choice of the operator affects the outcome of the joint inversion. The algorithm presented here offers the possibility to study such effects.

The performance of the 3-D structural joint inversion was assessed in a cross-well synthetic test case. The synthetic models, comprising two spherical anomalies in a layered background medium, were set up to mimic a scenario where on one hand the detection of the anomalies is a target, and on the other hand, the petrophysical distribution of the background medium is of interest. When joint instead of individual inversion was used, it was possible to obtain GPR-velocity and resistivity models that are closer to the true models. In particular, the ERT tomogram was enhanced using joint inversion. Moreover, the true petrophysical parameter distribution is better represented through the models from joint inversion. However, the full variability in the property fields could neither be obtained with joint nor with individual inversions. This is expected from regularized deterministic inversions and probabilistic inversions would be needed to retrieve the full variability. Furthermore, porosity distributions that are very similar to each other and show a good match with the true porosity distribution were calculated from the joint inversion results. For the individual inversions, the distributions are not as close to the true porosity distributions and differ from each other. In practice, a choice for one or the other porosity field needs to be made. While joint inversion produces nearly identical porosity fields, one needs to choose between two different fields for individual inversion.

Imposing prior geological information on the structural coupling operators yielded geologically meaningful joint inversion results. While dominant subsurface features are visible on both, individual and joint inversion results, the latter show more continuity of the dipping, layered structures, which are very typical for the survey area. Although we prefer the results from joint inversion, results from individually inverting the data should also be considered to asses the structural similarity and ensure that the models are not overly different, which might indicate that they are structurally not similar (Linde et al. 2008).

Appropriate weighting of the contributions associated with each method in the objective function is an issue with every joint inverse problem (Doetsch et al. 2010b; Lelièvre et al. 2012; Moorkamp et al. 2016; Steklova & Haber 2017). Our weighting approach is to use the data fit measure χ² to give more weight to the part in the inversion scheme (eq. 5) that has the higher data misfit. The motivation behind this approach is to achieve convergence in data fit at a similar pace for both methods. In both joint inversion examples presented here, this weighting scheme ensured similar convergence rates for the individual data misfits (see Figs 4 and 8).

7 CONCLUSIONS

We developed a structural joint inversion algorithm that imposes structural similarity between property fields over predefined physical length scales. In contrast to cross-gradient operators used in other studies, our operators are not only defined in the immediate neighbourhood of a cell but include a larger footprint that is defined via the correlation lengths of a geostatistical model. In particular for irregular meshes, where adjacent cells may vary in size, the larger footprint stabilizes the calculation of the cross-gradients. In a 3-D synthetic cross-well study, our joint inversion algorithm performed superior compared to individual inversions in discriminating anomalous bodies from a background model and recovering the petrophysical distribution of the GPR slowness and resistivity background models. For the joint inversion of field resistivity and traveltime data, prior geological knowledge was used to define the structural coupling and regularization operators. The resulting models are in accordance with the regional geology and known tectonic features are clearly imaged. Our methodology is applicable to structural joint inversion on any scale and will be especially beneficial for jointly inverting data from geophysical methods with strongly different subsurface resolution.

ACKNOWLEDGEMENTS

This work was supported by ETH Research Grant ETH-01-14-2. We thank the sponsors of the CARNEVAL consortium (Nagra, OMV and Schlumberger) for their support related to the acquisition of the field data sets. We thank Niklas Linde for fruitful discussions. Furthermore, we are grateful to the Swiss Institute for Speleology and Karst Studies (SISKA), in particular Urs Eichenberger, for sharing their expertise.