GEOPHYSICS, VOL. 79, NO. 6 (NOVEMBER-DECEMBER 2014); P. A45–A50, 5 FIGS.

10.1190/GEO2014-0220.1
Downloaded 10/17/14 to 129.11.21.2. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Inversion of gravity data with isostatic constraints

Ahmed Salem1, Chris Green2, Matthew Stewart3, and Davide De Lerma3

mentary basin. This enables the gravity anomaly to be modeled,

ABSTRACT especially if control data are available. Fairhead and Green
(1989) use a simple geometric model to represent the crustal thin-
We have developed a simple iterative gravity-inversion ning and model the residual gravity as being due to sedimentary rift
approach to map the basement and Moho surfaces of a rift basins. Jorgensen and Bosworth (1989) separate the gravity field
basin simultaneously. Gravity anomalies in rift basins com- using a best-fit, low-order polynomial to address a similar interpre-
monly consist of interfering broad, positive crustal-thinning tation. Salem et al. (2013) separate gravity data across the Red Sea
anomalies and narrow, negative sedimentary-basin anoma- into sediment and Moho effects based upon wavelength filtering. In
lies. In our model, we assumed that the Moho and basement each of these cases, the residual gravity can be inverted for base-
surfaces are in Airy isostatic equilibrium. An initial plane- ment depth using a method that inverts for the depth of a single
layered model was iterated to fit the gravity data. We applied interface, e.g., Cordell and Henderson (1968) in the space domain
the process to a model in which the inverted basement and or Oldenburg (1974) using fast Fourier transforms.
Moho surfaces matched the model surfaces well and to a In this paper, we present an alternative inversion method for grav-
gravity profile across the Kosti Basin in Sudan. ity data over simple rift models incorporating constant-density con-
trasts at the basement and at the Moho. We set a simple Airy-type
isostatic relationship as a constraint on our model, such that the
shape of the Moho and the base of the sediments are related through
INTRODUCTION the Archimedes principle. The method iterates from a plane-layer
starting model to a model that fits the data. We test our method on
Gravity anomalies from density variations within the crust are
theoretical gravity data over a basin model and field gravity data
often complicated in their shape and reduced in their amplitude
across a small rift basin.
by the effect of isostacy. A particular example is the case of
continental rift basins in which the thinning of the crust under
the rift causes a positive gravity anomaly, which tends to counteract
METHOD
the negative gravity anomaly of the sedimentary basin. Similar ef-
fects occur over and at the edges of ocean basins. Lateral strength in The basic isostatic model, sometimes known as the “hidden-
the crust gives rise to flexural isostacy, in which isostatic compen- layer” Airy model (e.g., Karner and Watts, 1982), which links
sation extends over a considerable distance from the near-surface the depth to basement hs and depth to Moho hm , is illustrated in
mass imbalance. However, in many stretched crust locations, it Figure 1 as
is found that compensation is relatively local and Airy isostacy
is a reasonable model, e.g., studies of the West African margin Δρs
hm ¼ hc þ hs ; (1)
by Watts and Stewart (1998). Δρm
Even for an Airy isostatic model, the two gravity anomalies do
not completely cancel as the deeper density contrast of the crustal where Δρs is the (negative) density contrast of the sediments rel-
thinning gives rise to a positive gravity anomaly, which is broader ative to the basement, Δρm is the (positive) density contrast at
than the negative gravity anomaly produced by the shallow sedi- the Moho (upper mantle relative to lower crust), and hc is the Moho

Manuscript received by the Editor 13 May 2014; revised manuscript received 16 July 2014; published online 1 October 2014.
1
Formerly Nuclear Materials Authority, Cairo, Egypt; presently GETECH, Kitson House, Elmete Hall, Elmete Lane, Leeds, LS8 2LJ, UK. E-mail:
[email protected]; [email protected].
2
Formerly University of Leeds, Leeds, UK; presently GETECH, Kitson House, Elmete Hall, Elmete Lane, Leeds, LS8 2LJ, UK. E-mail: chris.green@getech
.com.
3
GETECH, Leeds, UK. E-mail: [email protected]; [email protected].
© 2014 Society of Exploration Geophysicists. All rights reserved.

A45
 A46 Salem et al.

depth for basement depth of zero; all depths are positive downward. and gobs is below a predefined tolerance. The tolerance is defined on
The Bouguer gravity anomaly over the basic Airy model can be a case-by-case basis.
written as The convergence of the iteration defined by equation 3 is found to
depend on the iteration factor α. For the model example, α ¼ 1.0
gobs ¼ gm þ gs þ gc ; (2) iterates quickly to fit the data, with slightly higher iteration factors
reaching the solution more quickly. However, for the real-data ex-
where gm is the gravity effect due to the density contrast of the upper
Downloaded 10/17/14 to 129.11.21.2. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

ample, an iteration with α ¼ 1.0 does not converge to a good fit to

mantle, gs is the gravity effect due to the density contrast of the the data. In general, a lower value of α is more likely to give con-
sediment, and gc is a constant offset, possibly representing deeper vergence, albeit after a higher number of iterations.
density variation. The approach used in this study is similar to that For a realistic example, there is likely to be some regional com-
of Cordell and Henderson (1968), in which each gravity value has a ponent of the gravity field (represented by gc in equation 2), which
depth point in the model directly under it and the model is iterated must be removed before it is possible to invert for the basement and
until the gravity calculated from the whole model fits the data. In Moho surfaces. This regional component (the simplest case being a
this case, we consider a 2D model that continues to infinity in the uniform offset of the gravity data across the survey area) will pre-
perpendicular direction; each gravity-anomaly value lies directly vent the gravity model from converging to fit the gravity data pre-
above a model vertex in the basement surface and the Moho surface. cisely, so we remove it by adjusting the data profile during the
We iterate the model at each point independently based on the dif- iteration process. The regional adjustment can be approximated
ference between the observed gravity gobs and the gravity gcalc cal- based upon control data, but consecutive minor adjustments will
culated from the last result from the iterative inversion process. We still be required until the gravity data and the model fit well, because
use the “gpoly” algorithm (Blakely, 1995) for the forward calcula- the adjustment will be based on infinite slabs, whereas the models
tion of the gravity effect of the updated model. are 2D; interfering basement and Moho anomalies will also com-
The Bouguer slab model is used to inform the iteration, but be- plicate the situation. An alternative approach is to add a second iter-
cause the basement and Moho surfaces are linked by equation 1, we ation to the process in which the whole gravity observation profile is
cannot consider both these interfaces in the iteration because they adjusted by gadj to transform the previous gravity data values
would always balance. Instead, we focus on the upper (basement) ðgobs Þold to the new values ðgobs Þnew :
surface because it will have a bigger local effect at the gravity-ob-
servation point than the Moho surface. The iteration is based on the ðgobs Þnew ¼ ðgobs Þold þ gadj ; (4)
increase or decrease in the sedimentary thickness that would give a
gravity anomaly equal to the difference between the observed and where
previously calculated gravity anomaly. The basement depth is iter-
ated at each stage by 1X n
gadj ¼ ðg − gobs Þ (5)
n i¼1 calc
g − gcalc
ðhs Þnew ¼ ðhs Þold þ α obs ; (3)
2πGΔρs and n is the number of gravity observations. This gravity adjustment
is applied to the measured gravity data after each iteration of equa-
where ðhs Þnew and ðhs Þold are the depths to basement after and before
tion 3 to remove the regional component that cannot be removed
the iteration, respectively, G is the gravitational constant, and α is
based on the limited control points. This dual iteration is found
the iteration factor that controls the step size of the iteration process.
to converge rapidly to a stable solution.
The depth to Moho is then updated according to equation 1.
Once each point in the model has been iterated, the gravity effect
of the new model is calculated. A further iteration is defined based
MODELED-DATA EXAMPLE
upon the difference between the new gcalc and gobs . The process is
repeated until the root-mean-square (rms) difference between gcalc We created a gravity data set over a simple three-layer rift model
(Figure 2a) with basement, and Moho surfaces related by equation 1.
We used a density contrast of −400 kg∕m3 between sediments and
basement and a density contrast of 500 kg∕m3 between the upper
mantle and the crust. The gravity data were calculated along a 200-
km profile with a sample spacing of 1 km located above the model
nodal points. We added a value of 10 mGal to the calculated
anomaly, simulating a constant offset in the gravity anomaly. We
then contaminated the shifted gravity anomaly data with random
noise with zero mean and a standard deviation of 0.5 mGal.
Figure 2b is the resultant Bouguer anomaly profile that shows a
negative anomaly from the sedimentary basin superimposed onto a
positive anomaly from the crustal thinning; even though the positive
Figure 1. Physical constraints imposed by the Airy isostatic model, anomaly is far from complete in this window, we attempt to invert
showing the relationship between basement depth and Moho depth. for both the surfaces. We started with a simple plane-layer model
The Δρs value is the (negative) density contrast of the sediments
relative to the basement, the Δρm value is the (positive) density con- (Figure 2c), which necessarily fits equation 1. The model is consid-
trast at the Moho (upper mantle relative to lower crust), and the hc ered to extend horizontally to infinity at the ends of the model, and
value is the Moho depth for basement depth of zero. hence, the gravity anomaly calculated from this starting model is
 Inversion of gravity data A47

constant. We used the model values for Δρs, Δρm , and hc in the below the tolerance of 0.2 mGal, after 11 iterations in this case
inversion. We also used the depth to basement at a horizontal loca- (Figure 2d). Figure 2e shows the value of the adjustment in equa-
tion of 50 km as a depth control to remove an initial estimate of the tion 5 at each iteration. The iterative gravity adjustment started with
offset from the gravity data (12.9 mGal). The inversion process 3.0 mGal and reduced to −0.02 mGal after 11 iterations with a final
iterates from the initial model using equation 3 (iteration factor estimate of the offset from the gravity data of 9.7 mGal. Figure 2f
of 1.0) and equation 4 until the rms difference between the theoreti- shows that the final result for the test model fits very closely with the
Downloaded 10/17/14 to 129.11.21.2. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

cal gravity anomaly data and the calculated gravity anomaly of the basement and Moho surfaces, and Figure 2g shows the close fit of
updated model has been reduced from approximately 13 mGal to the calculated and original model gravity anomalies.

Figure 2. Inversion process for a rift basin model. (a) Three-layer model, (b) gravity anomaly from the model (every fourth point marked with a
symbol), (c) initial model for inversion, (d) the rms mismatch between the gravity from the model and the gravity from the result after each
iteration, (e) value of gravity adjustment after each iteration, (f) final inverted result of sediment/basement and crust/upper mantle interfaces
(circles) compared with the model (lines), and (g) gravity anomaly for the final result (lines) compared with the model (circles).
 A48 Salem et al.

FIELD-DATA EXAMPLE A starting model was created with densities of 2350 kg∕m3 for
the sediments, 2700 kg∕m3 for the crustal rocks, and 3300 kg∕m3
We demonstrate the practical utility of the approach using gravity for the upper mantle and hc ¼ 30 km — all as used by Jorgensen
data over the Kosti rift Basin in Sudan (Figure 3a). The Kosti Basin and Bosworth (1989). To estimate an initial regional offset from the
is located just to the west of the Blue Nile River and is one of the observed gravity data, we used a depth-control point at the shallow-
smallest of several northwest–southeast-trending Cretaceous rift ba- basement point closest to the basin (horizontal location 13 km). This
Downloaded 10/17/14 to 129.11.21.2. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

sins thought to have been formed as a result of lateral movement control point provided an initial estimate of the offset of approxi-
along the Central African Shear Zone (Browne et al., 1985; Jorgen- mately −23.3 mGal, which was then removed from the observed
sen and Bosworth, 1989; El Tahir et al., 2013). The basin is largely gravity data. To estimate appropriate densities, we used another
covered by surficial deposits of unconsolidated Cenozoic sands, control point (3.7 km depth at 24 km along profile) based on
gravels, shales, and clays with some Cretaceous clastic sediments the seismic interpretation of Jorgensen and Bosworth (1989).
and Precambrian basement outcropping to the northeast. The inversion process was iterated (with iteration factor 0.5)
A 65-km southwest–northeast gravity profile across the basin reducing the rms error between the observed and calculated gravity
was used to test the inversion method. This profile follows part to 0.17 mGal after 40 iterations (Figure 4b). The iterative gravity
of a longer profile interpreted by Jorgensen and Bosworth adjustment (Figure 4c) started with 16.2 mGal and reduced to
(1989). The purpose of using a limited profile length was to test −0.04 mGal after 40 iterations with a final estimate of the offset
the algorithm in a more challenging and realistic situation, in which from the gravity data of −38.9 mGal. The resulting model (Fig-
data gaps and interfering anomalies often make it impossible to in- ure 4a) fits the data well (Figure 4d) but has some large, short-wave-
vert the complete anomaly. The Bouguer anomaly across the Kosti length variations in the depth to the basement not seen in the model
Basin (Figure 3b) shows the classic rift basin anomaly pattern of a of Jorgensen and Bosworth (1989). The depth-to-basement result at
central low from the thick sediments superimposed on a broader the location of the seismic depth-control point is overestimated by
high from crustal thinning without apparent interference from other 1.9 km (Figure 4e); this discrepancy is seen to increase beyond iter-
gravity anomalies; as such, it should provide a good example for the ation five, which causes some concern about the inversion.
inversion method. The iterative inversion was rerun using a slightly different density
scheme: 2350 kg∕m3 for the sediments, 2800 kg∕m3 for the crustal
rocks, and 3300 kg∕m3 for the upper mantle. The rms error between
the observed and calculated gravity reduced to 0.16 mGal after 17
iterations (Figure 5b). The iterative gravity adjustment (Figure 5c)
started with 17.7 mGal and reduced to 0.02 mGal after 17 iterations
with a final estimate of the offset from the gravity data of
−39.0 mGal. This indicates that estimation of the gravity offset
is not affected by using a different density scheme.
The resulting model (Figure 5a) also fits the data well (Figure 5d).
The depth-to-basement results with the new density fit the seismic
control better (Figure 5e) and have less extreme depth variations.
The Moho depth model still has much greater variation than the
model of Jorgensen and Bosworth (1989); their Moho model, how-
ever, is very likely to be smooth because the regional gravity from
which it was modeled contains only very long wavelengths, having
been generated from low-order polynomials.

DISCUSSION
The application of the iterative inversion approach to generate
basement and Moho surfaces has generally been successful. A
few observations can be made. It appears that the starting model
is not particularly important and that the method will iterate to a
robust model that fits the data well, even from a very simple starting
model with no structure. However, the densities used in the model
are important; it appears from the field-data example (Figures 4 and
5) that changing the densities in the model slightly can generate a
much more (or less) realistic depth model. On the other hand, it
appears that the method could be diagnostic of whether the densities
used are correct, such that testing which densities give best conver-
gence could form part of a standard workflow. Different densities
can be assessed based on how rapidly the inversion reaches the re-
quired rms misfit and whether the control data misfit is large or
Figure 3. (a) Location of gravity profile A-A′ across the Kosti rift
Basin in Sudan. Basin outlines are modified from El Tahir et al. unstable. Better control data should help to identify correct densities
(2013). (b) Gravity anomaly along profile A-A′ (modified from Jor- earlier in the process; specific density information could be used, or
gensen and Bosworth, 1989). appropriate densities could be determined from depth-control points
 Inversion of gravity data A49

combined with gravity values. The gravity anomaly offset is also sion approach are inherently better or worse than smooth Moho sur-
important and significant offsets in the data can limit the conver- faces, such as those interpreted by Jorgensen and Bosworth (1989).
gence to a robust solution. A method has been developed to iter-
atively solve for the gravity offset, but this is based on the
CONCLUSION
assumption of a simple shift in the gravity data. If the regional grav-
ity field is more complicated, then removing the very longest wave- From tests on a modeled gravity data set and a field data set of a
Downloaded 10/17/14 to 129.11.21.2. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

lengths prior to the inversion might be required to generate good simple rift basin anomaly, we find that it is possible to invert for
results. basement and Moho surfaces together, provided the two are strictly
Applying the method in 3D is an obvious next step. Extension to linked with an equation based on the Airy isostatic model. This is
3D is conceptually straightforward. Some complications are likely despite the fact that the two gravity effects tend to counteract each
to occur in the details of how the regional field is defined or in how other. The method appears to consistently converge to a result that
the model is extended at the edge. The Moho is normally deep, and fits the gravity data very well and also replicates test models to great
gravity anomalies caused by it are generally long wavelength, mak- accuracy. This method appears to offer an alternative to approaches
ing it difficult to resolve its shape in detail from gravity data. In this that split the data into regional and residual fields using a math-
study, we impose a strict adherence to Airy isostacy, but unless we ematical approach — either wavelength filtering or polynomial
have access to other information, it will be difficult to decide separation. The new method does not make a mathematical
whether the rough Moho surfaces derived from this iterative inver- assumption but relies on a reasonable physical model of isostacy;

Figure 4. Inversion results for the Kosti Basin. (a) Basement and Moho surfaces. The dotted lines are the interpretation of Jorgensen and
Bosworth (1989). (b) The rms mismatch between the gravity from the model and the observed gravity after each iteration, (c) value of gravity
adjustment after each iteration, (d) gravity anomaly for the resultant model (lines) compared with the observed gravity data (circles), and
(e) absolute depth error at the location of the seismic depth-control point.
 A50 Salem et al.
Downloaded 10/17/14 to 129.11.21.2. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Figure 5. Same as Figure 4 but with a revised crustal density.

as such, it has the promise to produce more realistic models that Cordell, L., and R. Henderson, 1968, Iterative three‐dimensional solution of
gravity anomaly data using a digital computer: Geophysics, 33, 596–601,
inherently fit with physical principles. As yet, the method has only doi: 10.1190/1.1439955.
been tested in rather simple situations; application to more complex El Tahir, N., A. Nyblade, J. Julià, and R. Durrheim, 2013, Crustal structure
geologic problems will require a greater quantity of control data of the Khartoum Basin, Sudan: Tectonophysics, 593, 151–160, doi: 10
.1016/j.tecto.2013.02.032.
and/or a more complete understanding of the relationship between Fairhead, J. D., and C. M. Green, 1989, Controls on rifting in Africa and the
the different parts of the model. regional tectonic model for the Nigeria and East Niger rift basins: Journal
of African Earth Sciences (and the Middle East), 8, 231–249, doi: 10
.1016/S0899-5362(89)80027-2.
ACKNOWLEDGMENTS Jorgensen, G. J., and W. Bosworth, 1989, Gravity modeling in the Central
African Rift System, Sudan: Rift geometries and tectonic significance:
Journal of African Earth Sciences (and the Middle East), 8, 283–306,
We greatly appreciate the constructive and thoughtful comments doi: 10.1016/S0899-5362(89)80029-6.
of J. Rickett and three anonymous reviewers. We thank R. Blakely, Karner, G. D., and A. B. Watts, 1982, On isostacy at Atlantic-type
S. Campbell, and S. Mazur for their help in discussion and com- continental margins: Journal of Geophysical Research, 87, 2923–2948,
doi: 10.1029/JB087iB04p02923.
ments on this work. Oldenburg, D. W., 1974, The inversion and interpretation of gravity anoma-
lies: Geophysics, 39, 526–536, doi: 10.1190/1.1440444.
REFERENCES Salem, A., C. Green, S. Campbell, J. Fairhead, L. Cascone, and L. Moor-
head, 2013, Moho depth and sediment thickness estimation beneath the
Blakely, R. J., 1995, Potential theory in gravity and magnetic applications: Red Sea derived from satellite and terrestrial gravity data: Geophysics, 78,
Cambridge University Press. no. 5, G89–G101, doi: 10.1190/geo2012-0150.1.
Browne, S. E., J. D. Fairhead, and I. I. Mohamed, 1985, Gravity study of the Watts, A. B., and J. Stewart, 1998, Gravity anomalies and segmentation of
White Nile Rift, Sudan, and its regional tectonic setting: Tectonophysics, the continental margin offshore West Africa: Earth and Planetary Science
113, 123–137, doi: 10.1016/0040-1951(85)90113-1. Letters, 156, 239–252, doi: 10.1016/s0012-821x(98)00018-1.