of 20
Geophys. J. Int.
(2020)
223,
1899–1918
doi: 10.1093/gji/ggaa425
Advance Access publication 2020 September 08
GJI General Geophysical Methods
A Bayesian 3-D linear gravity inversion for complex density
distributions: application to the Puysegur subduction system
Erin Hightower
,
1
Michael Gurnis
1
and Harm Van Avendonk
2
1
Seismological Laboratory, California Institute of Technology,
1200
E California Blvd, Pasadena, CA
91125
, USA. E-mail:
ehightow@caltech.edu
2
Institute for Geophysics, University of Texas,
10100
Burnet Road, Austin, TX
78758
,USA
Accepted 2020 September 4. Received 2020 September 4; in original form 2020 May 30
SUMMARY
We have developed a linear 3-D gravity inversion method capable of modelling complex
geological regions such as subduction margins. Our procedure inverts satellite gravity to
determine the best-fitting differential densities of spatially discretized subsurface prisms in
a least-squares sense. We use a Bayesian approach to incorporate both data error and prior
constraints based on seismic reflection and refraction data. Based on these data, Gaussian
priors are applied to the appropriate model parameters as absolute equality constraints. To
stabilize the inversion and provide relative equality constraints on the parameters, we utilize a
combination of first and second order Tikhonov regularization, which enforces smoothness in
the horizontal direction between seismically constrained regions, while allowing for sharper
contacts in the vertical. We apply this method to the nascent Puysegur Trench, south of New
Zealand, where oceanic lithosphere of the Australian Plate has underthrust Puysegur Ridge and
Solander Basin on the Pacific Plate since the Miocene. These models provide insight into the
density contrasts, Moho depth, and crustal thickness in the region. The final model has a mean
standard deviation on the model parameters of about 17 kg m
–3
, and a mean absolute error on
the predicted gravity of about 3.9 mGal, demonstrating the success of this method for even
complex density distributions like those present at subduction zones. The posterior density
distribution versus seismic velocity is diagnostic of compositional and structural changes and
shows a thin sliver of oceanic crust emplaced between the nascent thrust and the strike slip
Puysegur Fault. However, the northern end of the Puysegur Ridge, at the Snares Zone, is
predominantly buoyant continental crust, despite its subsidence with respect to the rest of the
ridge. These features highlight the mechanical changes unfolding during subduction initiation.
Key words:
Gravity anomalies and Earth structure; New Zealand; Inverse theory; Statistical
methods; Subduction zone processes.
1 INTRODUCTION
Inverse methods have become increasingly popular for addressing
a number of problems in earth science, particularly for subsurface
mapping. Gravity inversion, for determining either the densities or
depths of bodies of known density in the Earth, has been an es-
tablished method of mapping the Earth’s heterogeneities for some
time, though often with emphasis on the non-linear approach. In
non-linear gravity inversion, the densities and density contrasts of
the subsurface bodies are assumed to be known and one solves for
the geometry of the source, usually in terms of depth to a partic-
ular interface. These inversions include either methods operating
in the spatial domain (Medeiros & Silva
1996
; Prutkin & Casten
2009
;Camacho
et al.
2011
) or those operating in the wavenumber
domain (Parker
1972
;Oldenburg
1974
;Parker
1995
; Chappell &
Kusznir
2008
; Cowie & Kusznir
2012
;Bai
et al.
2014
). However,
despite the Fourier method being one of the classical approaches to
gravity inversion, wavenumber methods are often less effective in
recovering a fully 3-D solution with multiple sources and complex
geometry (Bear
et al.
1995
;Geng
et al.
2019
).
With the linear method, the unknowns are the densities of a dis-
cretized array of subsurface rectangular prisms and iteration is not
required in order to reach model convergence, except in the case
of testing variations in model regularization or other constraints.
Solving for the 3-D density distribution also indirectly solves for
the depth to key interfaces, such as the Moho, because we can inter-
pret such boundaries from sharp transitions in density. While linear
gravity inversion is an established method (Bear
et al.
1995
;Li&
Oldenburg
1998
; Silva
et al.
2001
; Silva Dias
et al.
2009
; Barnoud
et al.
2016
; Welford
et al.
2018
;Geng
et al.
2019
), many of the
studies using it only do so for relatively simple geological geome-
tries, such as a single sedimentary basin, mafic intrusion or volcanic
C

The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. All rights reserved. For
permissions, please e-mail:
journals.permissions@oup.com
1899
Downloaded from https://academic.oup.com/gji/article/223/3/1899/5902856 by California Institute of Technology user on 03 November 2020
1900
E. Hightower, M. Gurnis and H. Van Avendonk
feature (Medeiros & Silva
1996
; Silva
et al.
2001
; Barnoud
et al.
2016
). Successful application of this method to crustal scale studies
and tectonic margins, with variable approaches to the implementa-
tion, also exist (Welford
et al.
2010
,
2018
;Geng
et al.
2019
), but
few have applied this method to subduction zones. Subduction mar-
gins posses a complicated juxtaposition of structure and rock types
and significant and sometimes sharp lateral variations in density, as
opposed to passive continental margins, which often exhibit a more
gradual change in structure and rock type that is more easily handled
by smoothed inversions. We construct a 3-D linear gravity inver-
sion for an active subduction zone, demonstrating the successful
application of this method to more complex density distributions
and bolstering the validity of this method and its use in tectonic
applications.
Inversion has the advantage of providing statistical feedback on
solution quality. Specifically, within a Bayesian framework, the ob-
jective is to determine the posterior distribution of a set of parame-
ters given prior distributions and likelihood functions that describe
how the data relate to those unknown parameters (Tarantola
2005
;
Aster
et al.
2013
; De La Varga & Wellmann
2016
; Wellmann
et al.
2018
). The Bayesian approach is particularly useful for geophysical
inverse problems, which are in principle ill-posed because they are
inherently non-unique. For example, gravity data cannot distinguish
between a narrow density anomaly at depth or a wider source near
the surface (Li & Oldenburg
1998
; Welford
et al.
2018
;Geng
et al.
2019
). Consequently, one must introduce constraints and
apriori
information in order to transform them into well-posed problems.
With the Bayesian formulation, we can account for both error in
the data and error in our prior information to reduce how that error
may be carried over into the final model, and we can quantify the
error on our final solution via the covariance and resolution op-
erators. The Bayesian approach we use here offers improvements
over traditional gravity inversion and modelling techniques, where
one usually removes the effect of the topography and the Moho and
analyses the residual. Such an approach requires assuming constant
layer densities when in fact those densities are often unknowns,
and it requires assuming a known Moho depth that has to manually
and iteratively be adjusted by the user. This makes it difficult to
fully incorporate lateral changes in density. The Bayesian approach
is more flexible and capable of handling complex 3-D geometries
because it allows us to constrain where the boundary is
most likely
to be based on seismic data and what the densities are
most likely
to
be, while allowing both to vary in accord with the gravity data, the
final boundary location being dependent on the differential density.
As such, we are able to draw conclusions about the 3-D density
distribution in a tectonic setting that would otherwise not be as
apparent with traditional forward or inverse gravity methods that
require harder constraints or restrictions.
There are a number of common constraints widely used in grav-
ity inversion, including inequality constraints, which specify the
lower and upper bounds of parameter estimates; absolute proxim-
ity constraints, which specify that model parameters must be close
to a specified value, based on geological information at particular
points; and relative equality constraints, which specify that the spa-
tial variation of the model parameter values must be smooth (Silva
et al.
2001
). Absolute proximity constraints are rarely used alone
because there is often not enough prior information available to con-
strain all model parameters. An exception would be the minimum
Euclidean norm, or similarly zeroth order Tikhonov regularization,
which requires all parameter estimates to be as close as possible to
null values. This type of regularization is biased towards a solution
with minimum density and tends to concentrate mass anomalies
toward the surface, which is not entirely physical or useful for our
interpretation of the subsurface.
Minimum structure inversion, however, is a commonly used
method (Last & Kubik
1983
; Li & Oldenburg
1998
; Farquharson
2008
), utilized by codes such as GRAV3D (Li & Oldenburg
1998
).
To overcome the inherent insensitivity of gravity to depth and thus
the tendency for the inversion to concentrate mass near the sur-
face, these methods often apply a depth weighting (Li & Oldenburg
1998
). Applying absolute proximity constraints and inequality con-
straints to specific regions of the model, however, overcomes the
need for a depth weighting (Welford
et al.
2018
;Geng
et al.
2019
).
While traditional inverse methods do allow for the adjustment of
smoothing parameters, bounds on densities and variable weighting,
they usually do so under hard constraints on predefined bound-
aries where the density is allowed to vary but the geometry of the
boundary remains unchanged (Li & Oldenburg
1998
; Welford
et al.
2018
). In contrast, the probabilistic approach offers more flexibil-
ity. Previous comparisons between such probabilistic methods and
approaches such as those used by GRAV3D (Welford
et al.
2018
)
highlight these distinctions as well, and we refer the reader to these
sources for a more in depth comparison. These comparisons show
that while each method has its advantages and disadvantages, a
probabilistic approach using sparse seismic Moho constraints may
not always lead to better results, particularly when there are sig-
nificant lateral variations in crustal thickness and composition, as
it tends to concentrate more unreasonable densities into different
parts of the model to compensate (Welford
et al.
2018
). In contrast
to previous applications of this probabilistic method (Barnoud
et al.
2016
; Welford
et al.
2018
;Geng
et al.
2019
), our approach directly
incorporates constraints on the interface depths
and
on composition
via the mapping of seismic velocities to density, not only at the lo-
cations of sparse depth to Moho constraints from seismic lines, but
interpolated throughout the model domain and weighted according
to the spatial extent of the prior data. We also propagate the error
on the seismic velocities into the density prior to ensure that the
densities obtained vary within a range that is consistent with the
error in the seismic velocities and that the seismic data does not too
strongly dominate the final model obtained by the inversion, such
that it remains predominantly resolved by the gravity. Moreover,
the Bayesian approach allows us to directly evaluate the error and
statistical validity of our results in a way that does not assume the
seismic data is the full truth.
Due to the non-uniqueness of gravity, however, even with abso-
lute proximity constraints, some sort of smoothing or stabilizing
functional is needed to produce a meaningful solution. This can
come in the form of relative equality constraints such as either first
or second order Tikhonov regularization, which spatially minimizes
the first or second derivative of the physical property, respectively.
Relative equality constraints by themselves have a tendency to pro-
duce a blurred but still valuable model of the density anomalies
(Portniaguine & Zhdanov
1999
; Silva
et al.
2001
). However, when
combined with absolute equality constraints, this inversion tech-
nique is often able to produce accurate representations of the source
geometry and density (Medeiros & Silva
1996
; Silva
et al.
2001
).
As such, our method uses a combination of absolute and relative
equality constraints in the form of Gaussian priors based on exist-
ing geophysical data and a combination of first and second order
Tikhonov regularization.
There is distinction in the literature between traditional regular-
ization methods and proper Bayesian approaches to inverse prob-
lems. Traditionally, regularization modifies the function relating
the data to the source of its signal, in an effort to eliminate the
Downloaded from https://academic.oup.com/gji/article/223/3/1899/5902856 by California Institute of Technology user on 03 November 2020
Bayesian 3-D linear gravity inversion
1901
unstable problem by replacing it with a similar stable one. This
often involves a penalty on the inversion that guarantees a unique
solution (Calvetti & Somersalo
2018
). The Bayesian approach, on
the other hand, by modelling the solution as a random variable,
allows one to use the exact function relating the data to its source
and offers the flexibility of obtaining multiple reasonable solutions,
as the final posterior model is in fact a probability distribution.
However, the non-uniqueness of gravity inversion in particular re-
quires some form of regularization. The regularization method that
best bridges the classical deterministic theory and the Bayesian ap-
proach is Tikhonov regularization because instead of modifying the
model function, it solves a minimization problem (Calvetti & Som-
ersalo
2018
). In that sense, Tikhonov regularization is essentially
a smoothness prior and can be implemented within a probabilistic
framework, allowing the inversion problem to remain Bayesian even
though it involves regularization.
2 METHODS
2.1 Calculation of forward gravity
We model the subsurface density and structure of a defined region
and its associated effect on the gravity by discretizing the subsur-
face into a finite number of rectangular blocks. The gravitational
attraction of each rectangular prism is calculated and then summed
to compute the gravity field. The gravitational attraction of a homo-
geneous right rectangular prism relative to an observation point on
the surface is given as in Turcotte & Schubert (
2014
)as

g
=
ρ
2

i
=
1
2

j
=
1
2

k
=
1
μ
ijk
[

z
k
arctan


x
i

y
i

z
k
R
ijk


x
i
ln
(
R
ijk
+

y
j
)

y
j
ln
(
R
ijk
+

x
i
)]
,
(1)
where

x
i
=
(
x
i
x
p
),

y
j
=
(
y
j
y
p
),

z
k
=
(
z
k
z
p
)and
μ
ijk
=
(
1)
i
(
1)
j
(
1)
k
.
is the density contrast of the prism, and

is
the universal gravitational constant.
x
p
,
y
p
and
z
p
are the coordinates
of the measurement point, and
x
i
,
y
j
and
z
k
are the coordinates of
the corners of the prism, where (
i
,
j
,
k
)
=
(1, 2).
R
ijk
is the distance
from the measurement point to a corner at
x
i
,
y
j
,
z
k
and is given by
R
ijk
=
(

x
2
i
+

y
2
j
+

z
2
k
)
1
/
2
.
The sum defines the geometry of the prism relative to the ob-
servation point and can be extended to the case of multiple prisms,
such that each prism in the domain has a single geometry coeffi-
cient for each gravity observation point. We invert gravity data at
N
observation points to obtain the best-fitting estimate of the densities
of
M
subsurface prisms, or
M
model parameters. Eq. (
1
) then results
in an
N
×
M
matrix
G
that describes the geometry of each prism
relative to each observation point times

. The gravity anomaly
at any observation point due to the combined attraction of all the
prisms is the product of this matrix and
, which is an
M
×
1
vector containing the differential density of each prism, expressed
as

g
=
G
.
(2)
2.2 Linear least-squares inversion
We adopt the method for linear least-squares inversion as given in
Aster
et al.
(
2013
) and Tarantola (
2005
). For
N
data points and
M
model parameters, where
g
i
(
m
) is the model prediction of the
i
th
datum (the

has been omitted for clarity), the least-squares misfit
is:
F
(
m
)
=
1
2
N

i
=
1
(
d
i
g
i
(
m
))
2
.
(3)
For a linear model such as that given in eq. (
2
), the model derivative
is independent of the model parameters, and our prediction can be
written directly as
Gm
. The Gauss–Newton solution of the model
parameters that minimizes the least-squares misfit in eq. (
3
) is thus:
m
=
(
G
T
G
)
1
G
T
(
d
)
.
(4)
The data
d
are the observed gravity anomaly values, and the model
parameters to be estimated are the differential densities of each
discretized block in the subsurface.
We accommodate data errors and prior constraints on the model
parameters in the inversion via a Bayesian approach. Bayes theorem
states that the probability of the model parameters, given the data, is
proportional to the product of (1) the probability of producing those
data with the model and (2) the probability of the model itself.
P
(
m
|
d
)
P
(
d
|
m
)
P
(
m
)
.
(5)
P
(
m
) is a prior that we use to restrict the model parameters to certain
values given our existing geological knowledge.
In including the data error in the least-squares solution, we make
the key simplifying assumption that the data are independent. In
the case of gravity, we are incorporating the relative attraction of
both adjacent and distal blocks of mass, and if the data are gridded
with some form of interpolation, then they are arguably not truly
independent. However, given the complexity of the problem and its
physical geometry, the interdependence of the data is difficult to
quantify and the simplifying assumption that the data are indepen-
dent is sufficient to perform the inversion. We assume each data
point can be represented by a Gaussian distribution with known
error such that we can define a new least-squares misfit:
F
(
m
)
=
1
2
N

i
=
1

d
i
g
i
(
m
)
σ
d
i

2
,
(6)
where we are now minimizing the difference between the known
and predicted gravity, given the error in the gravity data. From
Bayes Theorem, minimizing this new misfit
F
(
m
) is equivalent to
maximizing
P
(
m
|
d
). To incorporate the data error into the model
parameter solution, we define a diagonal and symmetric weight
matrix
C
d
with the data variance on the diagonal. The solution
becomes:
m
=
(
G
T
C
d
1
G
)
1
G
T
C
d
1
d
.
(7)
2.3 Tikhonov regularization
Linear least squares, even when using the generalized inverse or
the truncated generalized inverse to handle small singular values, is
often insufficient for many inverse problems due to non-uniqueness
and instability, especially for high-dimensional problems. Thus, a
form of regularization must be applied. We use a combination of first
and second order Tikhonov regularization, which stabilizes the in-
version and acts as a relative equality constraint on the values of the
model parameters. First order Tikhonov regularization minimizes
the square of the first spatial derivative of the model parameters
(i.e. the gradient), thus serving to flatten the solution. Second order
Tikhonov minimizes the square of the second spatial derivative of
the model parameters (i.e. the curvature) and hence smooths the
solution. Zeroth order Tikhonov, on the other hand, favors models
Downloaded from https://academic.oup.com/gji/article/223/3/1899/5902856 by California Institute of Technology user on 03 November 2020
1902
E. Hightower, M. Gurnis and H. Van Avendonk
that are small and is identical to applying a Gaussian prior with a
mean of zero and minimizing the square of the model parameter
values themselves.
As Tikhonov regularization is equivalent to applying a prior that
enforces either small values, flatness, or smoothness, we can derive
the regularized solution by adjusting the misfit equation to reflect
the additional minimization of the model parameters or their first or
second derivatives.
F
(
m
)
=
1
2
(
d
Gm
)
T
C
d
1
(
d
Gm
)
+
λ
2
(
Lm
)
T
(
Lm
)
,
(8)
L
is either the identity matrix, a first derivative finite difference
operator, or a second derivative finite difference operator for zeroth,
first, or second order Tikhonov regularization, respectively.
λ
is
a constant controlling the strength of the regularization. As the
misfit remains exactly quadratic with the addition of the Tikhonov
regularization term, the inverse problem remains linear, and the
weighted and regularized linear least-squares solution becomes
m
=
(
G
T
C
d
1
G
+
λ
2
L
T
L
)
1
G
T
C
d
1
d
.
(9)
For 3-D models, first and second order Tikhonov are implemented
using the sums of the finite-difference approximations to the first
or second derivatives in each direction, respectively. Because the
discretization of the grid can be different in the
x
,
y
and
z
directions,
we apply three different regularizations, with associated constants
α
for the
x
-,
β
for the
y
-and
ζ
for the z-direction. The derivation of
the Tikhonov regularization matrices is given in Appendix A. For
three-dimensions, the weighted Tikhonov regularized solution is
m
=
(
G
T
C
d
1
G
+
α
2
L
x
T
L
x
+
β
2
L
y
T
L
y
+
ζ
2
L
z
T
L
z
)
1
(
G
T
C
d
1
d
)
.
(10)
Without a flatness constraint in the far-field, abrupt density
changes at the edges of the model domain result in a classical
gravity edge effect. Consequently, to ensure mathematical stability,
we impose an infinite edge boundary condition, which allows the
gravity to smoothly continue off the edges of the model area. We
accomplish this condition by padding the domain with edge prisms
that are sufficiently long that they extend far beyond the edge of
the gravity grid (on the order of 1000 km for the regional problem
with which we test the method). We also enforce this condition
during the inversion by using first order Tikhonov regularization
with a strong regularization coefficient to minimize the difference
between the edge parameters and the adjacent values so that their
predicted densities are the same. Thus, we apply different orders
and strengths of Tikhonov regularization to the edges and the in-
terior of the model simultaneously. The interior of the model has
second order Tikhonov imposed in the horizontal directions to allow
for smooth continuity of density bodies in the subsurface, and first
order Tikhonov is applied in the vertical direction, as it is better
equipped to allow for sharp contacts between layers of rock, while
strong first order is applied on the boundary.
As before, this variable order Tikhonov regularization can be
achieved by redefining the misfit equation, where separate
L
ma-
trices apply different weights to different sets of model parameters
and different directions. The full Tikhonov regularized solution,
with boundary conditions applied, is
m
=
(
G
T
C
d
1
G
+
L
)
1
(
G
T
C
d
1
d
)
,
(11)
where
L
=
α
2
L
x
T
L
x
+
β
2
L
y
T
L
y
+
ζ
2
L
z
T
L
z
+
b
2
B
x
T
B
x
+
b
2
B
y
T
B
y
,
(12)
b
is the weight of the first order Tikhonov regularization applied
to the boundary condition.
B
x
and
B
y
are the regularization matri-
ces that apply the boundary conditions in the
x
and
y
directions,
respectively.
2.4 Priors
Meaningful solutions consistent with existing geological knowledge
are obtained by applying absolute equality constraints as Gaussian
priors. In this approach, each parameter is forced to be close to a
mean value but is allowed to vary within a specified range. Different
regions of the model domain can have different priors depending
on (1) what we suspect the densities of the rocks in those areas are
and (2) how confident we are in those values based on their location
relative to the other data we have. The prior on each parameter is
given by the Gaussian probability density function
P
(
m
k
)
=
1
σ
p
2
π
exp
(
1
2
σ
2
p
(
m
k
μ
p
)
2
)
,
(13)
where
m
k
is the estimated model parameter value,
μ
p
is the expected
value of that model parameter based on our prior information and
σ
p
is the standard deviation of the prior for that parameter.
As with the data error and Tikhonov regularization, we define a
new misfit by adding the exponential component of the Gaussian
prior to the existing misfit:
F
(
m
)
=
1
2
(
d
Gm
)
T
C
d
1
(
d
Gm
)
+
α
2
(
Lm
)
T
(
Lm
)
+
1
2
(
μ
p
m
)
T
C
p
1
(
μ
p
m
)
.
(14)
Defining the prior covariance operator
C
p
as an
M
×
M
diagonal
matrix with the variance of the prior on the diagonal, we arrive at
the final data weighted, Tikhonov regularized solution with prior
constraints
m
=
(
G
T
C
d
1
G
+
L
+
C
p
1
)
1
(
G
T
C
d
1
d
+
C
p
1
μ
p
)
,
(15)
where
L
is defined as in eq. (
12
). This is the final solution vector used
in our inversion. The row or column number of the elements along
the diagonal of
C
p
correspond to the index number of that model
parameter. Likewise,
μ
p
is an
M
×
1 vector for which each element
corresponds to the density of single prism. To apply different priors
to different model parameters, one need only use the coordinates
of the model parameter centroids within the desired region to find
the appropriate model parameter index and apply a value to that
element. If an element on the diagonal of
C
p
1
is zero, then no prior
is applied to that model parameter.
2.5 Quantifying error
A key advantage of a Bayesian approach is that it allows us to statis-
tically evaluate the solution, via the posterior covariance matrix
C
of
the model parameters and the resolution matrix
R
. The covariance
matrix is defined as the inverse of the Hessian:
C
=
(
G
T
C
d
1
G
+
L
+
C
p
1
)
1
.
(16)
Here the values of
m
estimated by the inversion are the centre-points
of the posterior Gaussian, and the diagonal values of the covariance
matrix
C
are their associated variances.
Downloaded from https://academic.oup.com/gji/article/223/3/1899/5902856 by California Institute of Technology user on 03 November 2020
Bayesian 3-D linear gravity inversion
1903
The resolution matrix is determined from the covariance matrix
(Tarantola
2005
):
R
=
I
CC
p
1
,
(17)
where
I
is the identity matrix. If the resolution matrix equals the
identity matrix, the model is fully resolved by the data. This par-
ticular formulation of the resolution operator primarily allows us
to distinguish between those parameters that are resolved by inver-
sion of the gravity data and those that are resolved by the prior.
Mathematically, this can be written as:
tr
(
I
)
=
tr
(
R
)
+
tr
(
CC
p
1
)
(18)
meaning the total number of model parameters is the sum of the
number of parameters resolved by the data and the number of pa-
rameters resolved by the prior information (Tarantola
2005
). Higher
resolution (values closer to 1) means those parameter values have
mostly been determined by the inversion—in other words, we have
learned something from the gravity that we did not know
apriori
.
On the other hand, low resolution (values closer to 0) means the val-
ues of those parameters are almost entirely attributed to the prior.
This is the case for regions of the model where the prior is very
strong, that is a very small prior variance.
Ultimately, solution quality is based on the mean absolute error of
the gravity and the mean standard deviation of the model parameters
as determined from the diagonal of the covariance matrix, as well as
visual inspection of the model to determine its geological reason-
ability. Even with relative and absolute equality constraints, gravity
inversion remains non-unique and there are a number of model so-
lutions that could fit the data. It is possible to obtain a solution that
minimizes the misfit as required but that still appears geologically
unreasonable and must be disregarded as the most likely poste-
rior distribution of densities. However, the regularization and priors
ensure enough stability in the model that with the appropriate regu-
larization parameters
α
,
β
and
ζ
, the model obtained is geologically
sound and in line with our standing geophysical knowledge.
3 SYNTHETIC TESTS
Estimating optimal regularization parameters is difficult for gravity
inversion. We use an iterative technique on a series of synthetic tests
to determine
α
and
ζ
values that produce (1) the best fit between
the predicted and observed gravity and (2) the most geologically
reasonable solution, which for the synthetic models, is a nearly
complete recovery of the known density distribution. We conduct
these synthetic tests on a simplified lower resolution model of a sub-
duction system. In all synthetic tests, we construct a density model,
compute the forward gravity as given by eq. (
1
) and add Gaussian
noise to the gravity using a similar standard deviation to that of the
data set we will later use (about 1.7 mGal). We invert this gravity
for a range of Tikhonov regularization parameters and orders, with
or without priors on specific sets of model parameters, while at-
tempting to recover the known density distribution and judging the
stability of the inversion.
The performance of the inversion when used with first and second
order Tikhonov is tested using a simplified synthetic 3-D model of
a subduction zone (depicted in representative cross-sections in the
bottom row of Figs
1
and
2
). We test various combinations of the
horizontal regularization coefficient
α
and the vertical regulariza-
tion coefficient
ζ
for the cases of only first order Tikhonov, only
second order Tikhonov, and a combination of second order in the
horizontal and first order in the vertical. For each of these cases, we
test four additional classes of constraints: no priors, priors enforced
only on parameters within the water column, priors enforced only
on parameters within the water and crustal layers, and priors on all
parameters, including the mantle. The prism size is about 17.5 km in
the
x
-direction, 22.5 km in the y-direction, and increases from about
206 to 2060 m from shallow to deeper depths in the
z
-direction. The
α
and
ζ
values tested range from 10
3
to 10
8
. There are a total
of 10, 648 model parameters and 22, 500 data points, yielding an
overdetermined system. The synthetic density model is constructed
with a seawater density of 1027 kg m
–3
, oceanic crustal density of
2900 kg m
–3
, sediment density of 2300 kg m
–3
, continental crustal
density of 2700 kg m
–3
and mantle density of 3300 kg m
–3
.Wede-
fine differential density,
, by subtracting the lateral average of
each layer from the true density of each prism in that layer. The prior
densities, when applied, match those differential densities. The stan-
dard deviation of the priors, when applied, are 5 kg m
–3
for seawater,
80 kg m
–3
for the sedimentary and crustal rocks and 100 kg m
–3
for
the mantle.
The results for these synthetic tests are summarized in Figs S1–
S4, which show gridded results for each combination of
α
and
ζ
in panels corresponding to the order(s) of Tikhonov regularization
used (panel rows) and the set of priors used (panel columns). Grey
regions demarcate
α
and
ζ
combinations where the regularization
strength is too low to produce stable results. The minimum of each
test for both the mean absolute error (MAE) on the gravity and the
MAE on the model parameters is plotted in each of these figures
as well. Fig. S1 depicts the MAE between the true gravity field
of the synthetic model and the gravity predicted by the recovered
density distribution. Changes in the gravity misfit are much more
dependent on the order of regularization than they are on the pres-
ence of a prior. For first order Tikhonov alone, the misfit increases
dramatically above
α
values of 10
4
because the model becomes too
flat to correctly reproduce the shorter wavelength variations in the
gravity field. For second order Tikhonov, stability is achieved at
ζ
values of 10
2
in cases with limited priors, above which the gravity
error remains reasonably low until
α
values of about 10
7
. For the
combination of first and second order Tikhonov, the error remains
reasonably low until an
α
value of 10
7
and between
ζ
values of
10
1
and 10
3
. The lowest error on the gravity amongst all the tests
is about 1.29 mGal, which is less than the noise level of 1.7 mGal,
and occurs for the case of first order Tikhonov with no priors for
α
=
10
0
and
ζ
=
10
1
. The lowest gravity error occurs for the case of no
priors because without priors the model is allowed to take whatever
shape it must, subject to the smoothness constraint, to fit the data,
again highlighting the inherent non-uniqueness of the gravity.
However, to achieve a geologically reasonable model, priors must
be applied. For the case of enforcing a prior on all parameters, the
minimum gravity error is still only 1.36 mGal, so the fit to the
gravity data is not compromised by adding priors, while the fit to
the true density model is dramatically improved. Fig. S2 illustrates
the MAE between the predicted model parameter values and the
true model parameter values of the known density model. For most
combinations of different regularization orders and priors, too small
of an
α
or
ζ
value and the regularization is not strong enough to
provide a smooth and continuous density distribution, yielding non-
physical fluctuations in the density values (Figs
1
and
2
, columns
1 and 2). Alternative cross sections with results using different
regularization strengths are shown in supplementary Figs S5 and
S6. For
α
values that are too large, the solution smooths over any
density variations almost entirely. For cases with no priors or limited
priors, the misfit decreases with increasing
ζ
, but for cases with more
priors, the misfit begins to increase again with larger
ζ
values after
Downloaded from https://academic.oup.com/gji/article/223/3/1899/5902856 by California Institute of Technology user on 03 November 2020