# Geophysics

- © 2008 Society of Exploration Geophysicists. All rights reserved.

## Abstract

The structural approach to joint inversion, entailing common boundaries or gradients, offers a flexible and effective way to invert diverse types of surface-based and/or crosshole geophysical data. The cross-gradients function has been introduced as a means to construct models in which spatial changes in two distinct physical-property models are parallel or antiparallel. Inversion methods that use such structural constraints also provide estimates of nonlinear and nonunique field-scale relationships between model parameters. Here, we jointly invert crosshole radar and seismic traveltimes for structurally similar models using an iterative nonlinear traveltime tomography algorithm. Application of the inversion scheme to synthetic data demonstrates that it better resolves lithologic boundaries than the individual inversions alone. Tests of the scheme on GPR and seismic data acquired within a shallow aquifer illustrate that the resultant models have improved correlations with flowmeter data in comparison with models based on individual inversions. The highest correlation with the flowmeter data is obtained when the joint inversion is combined with a stochastic regularization operator and the vertical integral scale is estimated from the flowmeter data. Point-spread functions show that the most significant resolution improvements offered by the joint inversion are in the horizontal direction.

## INTRODUCTION

To determine petrophysical properties, state variables, and structural boundaries, it might be necessary to combine information provided by models obtained from different geophysical data (e.g., Tronicke et al., 2004; Bedrosian et al., 2007). Interpretation of several individually inverted data sets can be illuminating, but the results usually are affected by the resolution limitations of each model.

For consistent interpretations of multiple geophysical models, it would be advantageous to have inversion tools that have similar formulations of the inverse problem regardless of the type of geophysical data being inverted. This would allow models to be coupled, as long as the data have comparable spatial support. By joint inversion, we refer to coupled models that are obtained by simultaneously minimizing a misfit function that includes the data misfit of each data type. Joint inversion can improve the resolution of each geophysical model and provide models that are consistent with one another and therefore easier to interpret (e.g., Gallardo and Meju, 2004).

Joint inversion is not yet a standard tool in geophysical applications, mainly because robust and well-established petrophysical models that can be used to couple the models are usually available only for certain geophysical parameters, such as compressional and shear waves speeds (Tryggvason et al., 2002). Furthermore, petrophysical models often apply only in restricted geologic settings (e.g., Sen et al., 1981; Marion et al., 1992). In addition, the parameters of petrophysical models seldom can be constrained adequately by individual field data sets, so that fairly strong assumptions are required to couple models based on their petrophysical properties.

To avoid introducing questionable petrophysical models, joint inversion methods have been developed for layered (1D) structures that are expected to have coincident layer boundaries and constant properties within each layer (e.g., Monteiro Santos et al., 2006). A natural extension to 2D and 3D applications has been to assume that the earth can be divided into subvolumes of uniform properties with geometries that are common for all physical properties in the inversion (e.g., Hyndman and Gorelick, 1996; Musil et al., 2003). Such approaches are certainly useful, but physical properties can vary gradually in space, and not all data necessarily are sensitive to the same changes in lithology and state variables. Furthermore, the zonations must be updated continuously, making the inversions computationally expensive.

In Occam's inversion (Constable et al., 1987), fine model discretizations are used. The inverse problem is regularized by minimizing, for example, model roughness with the constraint that the simulated model response is close to a given target data misfit. Haber and Oldenburg (1997) introduce a joint inversion scheme to find models that are structurally similar, in the sense that spatial changes in models occur at the same location. This scheme is applicable to overparameterized 2D and 3D models. In essence, it is based on minimizing the squared difference of a weighted Laplace operator of the two models.

Gallardo and Meju (2003) further developed the framework of the structural approach to joint inversion proposed by Haber and Oldenburg (1997) by defining the cross-gradients function *x*, *y*, and *z*, and × indicates the cross product. By forcing the discretized cross-gradients function to be close to zero at each location and in each direction during the inversion process, either the gradients of the two resulting models will be parallel or antiparallel to each other or one or both of the models does not change. The boundaries of the resulting models have the same orientation, thus facilitating geologic interpretations. An advantage compared with the method of Haber and Oldenburg (1997) is that constraints based on the cross-gradients function do not focus on the magnitudes of the changes, which are difficult to estimate a priori and might necessitate some tuning parameters but in a common direction. The validity of imposing the cross-gradients function to be zero is discussed by Linde et al. (2006a) in the context of electrical properties.

The cross-gradients function first was employed in 2D for joint inversion of surface-based electrical resistance tomography (ERT) and seismic refraction profiles collected over weathered granodioritic bedrock overlain by mudstone (Gallardo and Meju, 2003; 2004). The same seismic refraction data later were inverted jointly with controlled-source audio magnetotelluric data (Gallardo and Meju, 2003). Recently, Gallardo (2007) proposed an extension of the cross-gradients function that allows simultaneous inversion of more than two data types. Gallardo (2007) then used this scheme in two dimensions to jointly invert surface-based P-wave and S-wave traveltimes, ERT, and magnetic data.

Tryggvason and Linde (2006) presented the first 3D application of joint inversion based on the cross-gradients function. By jointly inverting P- and S-wave traveltimes in a synthetic local earthquake tomography experiment, they found that the resulting anomalous

Linde et al. (2006a) presented the only study so far in which geophysical models obtained from joint inversions using the cross-gradients function were evaluated against borehole data. However, their example provided by the unsaturated Sherwood Sandstone was limited by the 1D character of the geology. We acknowledge that the cross-gradients method needs to be applied to additional well-instrumented study sites before it is accepted widely.

The objective of our study was to test whether a joint inversion could improve the resolution of lithologic boundaries compared with individual inversions of the same data sets. Here, we present the results of jointly inverted crosshole radar and seismic traveltimes recorded at the South Oyster study site in Virginia (Chen et al., 2001; Hubbard et al., 2001), where the geology consists of saturated, unconsolidated sediments with 3D heterogeneity. To guide the choice of inversion parameters, we calculate trade-off curves between the weight given to the cross-gradients constraint in the objective function and the resulting structural similarity of the resulting model as well as between the weight given to the cross-gradients constraints and the weight given to the regularization operators. Finally, the improved resolution offered by the joint inversion approach is visualized with point-spread functions calculated in the central part of the interwell region.

## METHODS

### Inversion method

Our formulation of the joint inverse problem closely follows that developed by Linde et al. (2006a) for use with ERT and GPR data. In the algorithm, the first-arrival traveltimes are computed with the finite-difference (FD) algorithm *time3d* (Podvin and Lecomte, 1991; Tryggvason and Bergman, 2006). Ray tracing is performed by a posteriori back propagation perpendicular to the wavefronts from the receivers to the transmitters (Vidale, 1988).

Our objective function has three competing components: the data fit resulting from two piecewise constant individual slowness models *λ* that scales the predicted cross-gradients function for a proposed model. A high *λ* gives a large relative weight to structural similarity in comparison with the regularization and the data misfit terms. These objectives are formulated in one objective function that we seek to minimize, at each iteration, in a least-squares sense (see equations 8–9 in Linde et al., 2006a) with the iterative conjugate gradient algorithm LSQR (Paige and Saunders, 1982). Individual inversions are performed by considering only one data set and model at a time (see equations 1–2 in Linde et al., 2006a).

We regularize the inverse problem by seeking either a model that has small second derivatives (i.e., smoothness constraints) or is close to a given stationary exponential covariance function. Linde et al. (2006a) show how such stochastic regularization operators can be calculated efficiently. To include the structural constraints in our objective function, we discretize the cross-gradients function, equation 1, with a central finite-difference scheme, which provides better results than the forward-difference scheme used in previous work (e.g., Gallardo and Meju, 2004). The cross-gradients function for a candidate model is estimated using a first-order Taylor expansion around the cross-gradients function at the previous iteration (see equation 7 in Linde et al., 2006a).

### Inversion strategy

The goal of the inversion is to construct two models so that each model satisfies the data within predefined data errors (i.e., the weighted rms of each of the models is very close to 1, indicating that the data are fitted to their estimated error levels) for a given *λ*, with the smallest possible values of *λ* fixed during the inversion and choose it so that the objective function is dominated by the terms that force the linearized cross-gradients function to be close to zero (i.e., *λ* is large). In the section on joint inversions, we investigate how the choice of *λ* affects the resulting model.

The joint inverse problem is nonlinear, not only because the raypaths depend on the slowness structure but also because the cross-gradients function is nonlinear (see equation 1). The derivatives of the cross-gradients function are valid only in the vicinity of the previous model, and care must be taken that the model updates are made within the region where the linearized cross-gradients function is reasonably valid.

To avoid large model updates outside the region where the linearizations are valid, we keep

To visualize the results, we introduce a normalized cross-gradients function as
*y*-component

## SYNTHETIC EXAMPLE

A synthetic example is used to demonstrate the impact of the cross-gradients constraints on the solution of the inverse problem. Figure 1a shows a GPR velocity model that consists of three homogeneous rectangular zones of high

Our synthetic experiment involved seismic and GPR pulses transmitted at intervals of

The resulting radar and seismic velocity models of the individual inversions are shown in Figure 1b and e. The models are smeared, making it difficult to determine the actual geometry of the zones. The models also indicate that the locations of the highest velocities in the common high-velocity zone differ for the seismic and radar models. The radar and seismic velocity estimates are plotted for every collocated model pair between boreholes (Figure 2a). This scatter plot seems to indicate that the subsurface is composed of three anomalous zones.

The resulting radar and seismic velocity models of the joint inversion are shown in Figure 1c and f. The three zones are less smeared out, and the geometries of the upper and lower zones correspond well with the actual geometries. In addition, the common high-velocity zone in Figure 1c and f has fairly uniform velocities. The scatter plot in Figure 2b shows very little scatter in the direction normal to the main axes of the velocity gradients, and the anomalies therefore can be picked by eye. For additional examples and discussions illustrating this situation, see Gallardo and Meju (2004, 2007), Linde et al. (2006a), and Tryggvason and Linde (2006). All models have a weighted rms very close to one.

As expected, the models constructed from the individual inversions and the joint inversion all resolve the high-velocity zones better than the low-velocity zones. Figure 2b indicates that the joint inversion resolves the upper zone of high seismic velocity better than the corresponding zone of low radar velocity, and that it resolves the lower zone of high radar velocity better than the corresponding zone of low seismic velocity. However, the joint inversion improves the geometry of the low-velocity zones (e.g., compare the low-velocity zones in Figure 1b and c).

The velocity contrasts considered in the example above are larger than the ones encountered in most environmental applications. Figure 2c and d shows the scatter plots resulting from individual inversions and joint inversion for a model in which the velocity range of the radar model was

## FIELD EXAMPLE

We now consider radar and seismic traveltime data collected between wells S14 and M3 in the South Oyster focus area, Virginia (Hubbard et al., 2001). Radar data were collected using a PulseEKKO 100 system with 100-MHz nominal-frequency antennae and a transmitter and receiver spacing of

We used a cell-size discretization of

### Individual inversions

Individual inversions were carried out using two types of regularization. The first was an anisotropic roughness operator that penalizes spatial variations in the horizontal direction five times as much as in the vertical direction. The second was a stochastic regularization operator (Linde et al., 2006a) based on an exponential model with a vertical integral scale of

The individually inverted smoothness-constrained radar (Figure 3a) and seismic (Figure 3e) velocity models contain a common low-velocity zone in the upper

The individually inverted radar (Figure 3b) and seismic (Figure 3f) velocity models based on stochastic regularization include the same main zones as those based on smoothness-constrained inversions. However, the former models are more variable because stochastic regularization operators represent a combination of general smoothness constraints (i.e., providing smooth models) and damping constraints (i.e., providing uncorrelated small-scale variability around an a priori model) (Maurer et al., 1998). The short vertical integral scale used here makes the zones thinner, more variable, and more pronounced. Models based on stochastic regularization tend to split the upper low-velocity zones into two horizontally aligned zones. A similar behavior is found in the high-velocity zone at

### Joint inversions

A series of joint inversions with stochastic regularizations was carried out with different weights *λ* given to the cross-gradients constraints (see equations 8–9 in Linde et al., 2006a). Figure 4a indicates that the smallest value of the cross-gradients function is obtained for a *λ* of 10,000, which results in models in which the cross-gradients function is approximately 100 times smaller than for the individual inversions. Figure 4b shows that a *λ* in the range of 100 to 10,000 causes the trade-off parameter in the final iteration *λ* grows above 10,000, *λ* of 10,000, corresponding to models that are structurally the most similar and with spatial variabilities that are only slightly larger than those of the individual inversions.

The models obtained by jointly inverting radar (Figure 3c) and seismic (Figure 3g) traveltimes with anisotropic smoothness constraints contain the same main zones as in the individually inverted models. However, the resolution is higher because more small-scale variability is present and the boundaries between the different zones are sharper. The low-velocity zone in the upper

The associated scatter plot (Figure 3o) is less scattered and thereby easier to interpret than those of the individually inverted models (Figure 3m and n). As expected, a strong positive correlation exists between the radar and seismic velocities. However, there are also zones in which the relationship changes, thus demonstrating the flexibility of the structural approach to joint inversion in dealing with nonstationary apparent petrophysical relationships. The correlation coefficient is 0.88.

The models obtained by jointly inverting radar (Figure 3d) and seismic (Figure 3h) traveltimes with stochastic regularization operators are spatially the most variable, but they contain mostly the same features as observed in models based on anisotropic smoothness constraints. The associated scatter plot (Figure 3p) indicates a strong linear correlation between the radar and seismic velocities, with a correlation coefficient of 0.90.

The values of the resulting cross-gradients function are approximately 100 times smaller for the joint inversion using the smoothness constraints (Figure 3k) and the stochastic regularization (Figure 3l) compared with the corresponding individual inversions (Figure 3i and j). The zones of high and low cross-gradients values are also much smaller in the joint inversions than in the individual inversions.

Table 1 demonstrates that the rms data fits for all models are quite similar. The joint inversions decrease the spread in the scatter plots significantly (Figure 3o and p), but they do not markedly change the trends that appear in the scatter plots for the individual inversions (Figure 3m and n). This result indicates that the regularization operator also is very important, because it will influence any zonation or petrophysical interpretation of the models significantly. The visual similarity between the individual inversions and the joint inversions indicates that the joint inversion is not getting trapped in local minima.

### Image appraisal

To investigate the resolution characteristics of the resulting models, we calculate the point-spread function (PSF) for the final inversion models based on the stochastic regularizations. The PSF is a row of the resolution matrix (e.g., Alumbaugh and Newman, 2000). It can be interpreted as the spatial averaging filter that relates the true underlying model to the resulting model estimate at a specific location. The PSFs are calculated with the LSQR algorithm in a method analogous to the one presented by Alumbaugh and Newman (2000) for the conjugate gradient method and with the assumption that the cross-gradients function for the true model is zero. The PSFs calculated for the joint inversions are normalized first with regard to the mean values of the radar and seismic slownesses. The PSFs then are normalized with regard to the largest value, in accordance with Alumbaugh and Newman (2000).

Horizontal and vertical profiles through the normalized PSFs at a central location (

The radar model from the joint inversion is dependent at this location only on the surrounding radar model and not on the seismic model (Figure 5a and b). The seismic model is affected here primarily by the neighboring seismic model, but it also is affected strongly by the surrounding radar model. This observation can be explained by studying the models at this location. The radar model indicates a low-velocity zone with a significant gradient that influences the resulting seismic model through the cross-gradients constraints (see equation 7 in Linde et al., 2006a). The seismic model is fairly uniform at this location and therefore does not impose any significant restrictions on the radar model. The PSFs constructed for the joint inversion are thus strongly dependent on the neighboring values of the seismic and radar slowness models. This means that the PSFs are meaningful only if they are evaluated around the final model estimate.

### Comparison with flowmeter data

Hubbard et al. (2001) used tomograms from the Oyster site along with flowmeter data to create a hydraulic conductivity model. The flowmeter data collected in the boreholes were kriged first to provide a prior model of hydraulic conductivity. Relationships between geophysical model parameters and hydraulic conductivity were developed using collocated tomographic estimates and flowmeter data. These relationships had a linear correlation coefficient of 0.68 (between hydraulic conductivity and radar velocity) and 0.67 (between hydraulic conductivity and seismic velocity). A Bayesian model then was used to update this prior model based on the observed correlation between flowmeter data and collocated geophysical model parameters. Scheibe and Chien (2003) found that model predictions based on stochastic indicator simulations conditioned to these estimates were superior to those based on stochastic indicator simulations conditioned to the flowmeter data alone.

The usefulness of geophysical tomograms for constraining hydrologic models depends strongly on the intrinsic relationships between the geophysical properties and hydraulic conductivity (e.g., Linde et al., 2006b), here most likely through a common link with porosity (e.g., Carcione et al., 2007). If the models based on the joint inversions correlate better with the flowmeter data, it should be possible to improve the hydraulic conductivity model of Hubbard et al. (2001) and therefore also the solute transport predictions.

To evaluate whether the joint inversions are likely to provide better models than those supplied by the individual inversions, we have compared hydraulic conductivity estimates based on flowmeter data at borehole M3 (located near the right side of the tomographic models in Figure 3) with tomographic estimates located two model cells away from the borehole. The flowmeter data are shown in Figure 6a, and the collocated radar velocity and seismic velocity models are shown in Figure 6b and c. To facilitate comparison, we highlight three depth intervals in which the hydraulic conductivity data display a local minimum (L1-L3) and three in which they display a local maximum (H1-H3).

The radar velocity model based on individual inversion with stochastic regularization accurately indicates five zones, whereas the radar model based on anisotropic smoothness constraints accurately indicates only one zone. The seismic velocity model based on stochastic regularization accurately indicates five zones, whereas the seismic model based on anisotropic smoothness constraints does not indicate any of those zones. The correlation coefficients of the collocated models based on the separate inversion and the flowmeter data are higher for the stochastic regularization (0.72 for radar velocity and 0.60 for seismic velocity) than for the anisotropic smoothness constraints (0.63 for radar velocity and 0.49 for seismic velocity).

The models based on the joint inversion with stochastic regularization indicate all six zones. The radar model based on anisotropic regularization indicates three zones, and the seismic velocity model indicates four zones. The correlation coefficients of the collocated models based on joint inversion and flowmeter data are higher when using stochastic regularization (0.78 for radar velocity and 0.69 for seismic velocity) than when using anisotropic smoothness constraints (0.67 for radar velocity and 0.52 for seismic velocity).

These results, which are summarized in Table 2, indicate that joint inversion with stochastic regularization at this site offers higher resolution and accuracy in determining lithologic changes than do the other inversion approaches considered here. Stochastic regularization operators based on borehole data (flowmeter data) provided better models than those supplied by traditional anisotropic smoothness constraints. Indeed, changes associated with varying the regularization operator seem to improve the resulting models slightly more than the joint inversion. Individually, the radar and seismic traveltime data could not provide the necessary resolution to image the small-scale variability revealed by the flowmeter data. Regularization based on flowmeter data from the site and the assumption that the two models are structurally similar improves the resolution of the resulting models.

## DISCUSSION

The resulting inversion models are qualitatively similar, and they reach the target misfit for all inversion types considered. This indicates that the nonlinearity of the cross-gradients function does not cause the joint inversion process to be trapped in local minima for the examples considered here. When performing this type of joint inversion, we recommend first to perform individual inversions to assess whether the models appear to be structurally similar. In addition, the individual inversions help to define a target data misfit and to assure that

The best way to choose *λ* when performing the joint inversions is to perform a series of inversions with different choices of *λ* (see Figure 4). The resolution improvements offered by the joint inversion vary from case to case and throughout the model domain because they are dependent largely on the spatial variations of the physical parameters we wish to recover. To evaluate the resolution improvements, we recommend calculating PSFs at selected locations (see Figure 5), as suggested by Alumbaugh and Newman (2000). We expect that some of the oscillations in the horizontal profiles of the PSF (e.g., Figure 5c) are caused by the discretization of the cross-gradients constraints, which operates on neighboring model cells only. This problem could be avoided if the cross-gradients constraints were discretized on the same scale as the stochastic regularization operator.

The stochastic regularization operator defines a presumably known scale on which structure is resolved. Large integral scales will result in large-scale features with large amplitudes, whereas small integral scales will result in a heavily damped solution with much small-scale variability (see Maurer et al., 1998, for an illustration). This implies that the magnitudes exhibited in the scatter plots of radar and seismic velocities (e.g., Figure 3p) are dependent on the regularization operator. The synthetic example also shows that earth models with large variability (e.g., Figure 1a and d) will introduce bias in the resulting scatter plots (see Figure 2b) because the velocity of low-velocity zones will be overestimated. These effects must be considered when making petrophysical interferences from such scatter plots (e.g., Linde et al., 2006a).

The sharper boundaries and the more variable models offered by the joint inversion with stochastic regularization are likely to improve flow and transport predictions significantly if used to constrain hydrologic models compared with using individual inversions with smoothness constraints. In the future, we plan to develop new strategies to infer the underlying geometry and physical properties of the earth from geophysical models resulting from joint inversion. The resulting method will be used to parameterize a 3D flow and transport model at an active research site.

## CONCLUSIONS

We used the cross-gradients function to develop an iterative nonlinear traveltime tomography algorithm that jointly inverts crosshole seismic and GPR traveltimes. From a synthetic example, we conclude that for structurally similar models, the joint inversion improves the determination of lithologic boundaries compared with models based on individual inversions.

The inversion method was applied to crosshole GPR and seismic traveltime data collected in an unconsolidated, saturated environment. The radar and seismic velocity models based on individual inversions indicated a strong structural similarity. The magnitudes of the cross-gradients function for the models based on the joint inversion were approximately 1% of the individually inverted models. The jointly inverted models were not only structurally more similar, but they also displayed larger spatial variability for a given data fit.

We used flowmeter data collected in one of the boreholes to assess whether the stochastic regularization and the joint inversion could improve the resulting models as indicators of geologic variability. We identified three zones in which flowmeter data indicated local minima in hydraulic conductivity and three zones indicating local maxima. The joint inversion based on the stochastic regularization was the only one that accurately located all six zones. Because tomographic models can be used to constrain hydrologic flow and transport models, we conclude that if structural similarity can be established, then the structural approach to joint inversion and the use of borehole data to determine regularization operators potentially can improve such hydrologic models.

We recommend this approach to joint inversion (1) when models based on individual inversions are structurally similar and (2) when structural similarity of geophysical properties can be expected from petrophysical considerations.

## ACKNOWLEDGMENTS

We thank Alan Green and Stewart Greenhalgh for commenting on draft versions of this manuscript. We also thank associate editor Lanbo Liu and three anonymous reviewers for constructive reviews that helped to improve the paper. Funding for this study was provided partially by the Swedish Geological Survey under contract 60-1307/2004 and by the Environmental Remediation Science Program, Office of Biological and Environmental Research, U. S. Department of Energy Grant DE-AC03-76SF00098.

- Received September 23, 2007.
- Revision received March 19, 2008.