NSW Isostatic Gravity - Roberts Geospatial

Correcting the New South Wales Gravity Database: Bouguer, Terrain, Spherical Cap, and Isostasy Corrections

Sam Roberts, Roberts Geospatial Engineering, April 2024

Introduction

The Australian National Gravity Database (ANGD), clipped to New South Wales (NSW), comprises 144,945 gravity station observations (correct on 2 April 2024). In this unsponsored study, I processed a spatially thinned subset of this database comprising 26714 (18%) observations. I selected these observations by spatially searching for the closest observation to a 0.05-degree grid covering NSW.

A correction was computed for each station. This correction is the vertical gravitational component of the gravitational field of an idealised Earth model constructed for each station and centred on that station location. This gravitational response is subtracted from the Free Air Gravity, and an arbitrary DC shift of 24,580 milligals was added to bring the median of the corrected gravity to approximately zero.

An Earth Model for each station is generated in suitable projected (UTM) coordinates, using the UTM zone in which the station is located. The model, viewed top-down, is circular and extends to a radius of 819,200 metres around the station. It is based on a spherical cap approximation of the Earth, so it curves downwards spherically away from the station. Onto this spherical curve is superimposed terrain to create a model of the topography of the Earth.

The Earth Model is comprised of three bodies representing the Continental Crust, the Mantle, and the Ocean. It incorporates onshore terrain, offshore bathymetry, and isostasy depth variations at the crust-mantle boundary. It effectively combines the Bouguer Correction, the Terrain Correction, the Spherical Cap Correction, and the Isostatic Correction – all into one accurate gravity correction.

The illustration shows the onshore and offshore topography of Southeast Australia. The major features include the Great Dividing Range running down the East Coast, the abrupt termination of the continental shelf along the East Coast, and transition to deep water and oceanic crust.

Spherical Cap Bouguer + Terrain Correction

I have previously described software that I have developed to compute a Spherical Cap Complete Bouguer Correction for gravity survey data. You can find more information here:

https://robertsgeospatial.com.au/spherical-cap-complete-bouguer-correction-processing-methodology/
https://robertsgeospatial.com.au/spherical-cap-complete-bouguer-correction-traverse-examples/
https://robertsgeospatial.com.au/spherical-cap-complete-bouguer-correction-grid-examples/

You can download documents describing the technique here:

Documents

Originally, the software generated a model of the Earth at each survey location, based on a spherical cap approximation of the Earths surface, that incorporated onshore terrain. The model was truncated at the shoreline and did not incorporate any offshore islands or other terrain. The model extended downwards to the datum plane, usually defined as 0m RL. When the vertical gravity component of this model is computed, it replaces the Bouguer Correction, the Spherical Cap Correction, and the Terrain Correction. I call this the Spherical Cap Complete Bouguer Correction.

An important point to note is that the terrain correction is computed for a spherical Earth, not a flat Earth. This is different to all other implementations and improves the quality and correctness of the terrain corrections. In a flat Earth model, both mass above the measurement plane (a mountain) and missing mass below the measurement plane (a valley) reduce the measured gravity. In both cases the terrain correction is additive. In a spherical Earth model, the further you move away from the measurement location, the further the Earth plunges below the measurement plane. In the far zone, a mountain may lie in-part or wholly below the measurement plane and consequently increase the observed gravity, reversing the sign of the terrain correction. An additional effect is that the horizontal distance from the measurement location to mass in the far zone is less than in a flat Earth model, and the vertical distance is increased. Both factors impact the terrain correction.

Offshore Bathymetry

Where gravity surveys are conducted near the coastline, the gravitational effect of the ocean and of offshore bathymetry, which can dwarf onshore terrain, can be considerable and it is desirable to correct for this effect.

The software was extended to create two mutually exclusive models at each station. The first is the Continental Crust and the second is the Ocean. The top surface of the crust model is the air-rock boundary onshore and the water-rock boundary offshore. The bottom surface is at the datum plane. This was chosen to lie below the ocean floor. The top surface of the ocean model is sea level (0m RL) and the bottom surface, derived from bathymetry, is the water-rock boundary. The crust and ocean models share a common boundary (water-rock).

Isostasy

The software has now been extended to incorporate isostatic effects and is now comprised of three models at each station. To paraphrase Wikipedia, isostasy is a state of gravitational equilibrium between the Earth’s crust and mantle such that the crust “floats” on the mantle at an elevation that depends on its thickness and density. For more information see:

https://en.wikipedia.org/wiki/Isostasy

To incorporate isostasy in the Earth model, I vary the bottom surface of the crust model, which is mirrored by the top surface of the mantle model. The variations in the depth of the crust are computed using the Airy–Heiskanen model of isostasy. In this model the crust has a constant density, and the thickness of the crust varies as the terrain height varies. Given the density of the crust, the mantle, and the ocean, the thickening of the crust below mountain ranges and the thinning of the crust in ocean basins can be computed.

Constructing the Earth Model

The Earth model consists of three bodies. Each body is a closed spherical cap with top and bottom surface variations. The surface of the body is a perfect triangular facet mesh with no holes or cracks. Where the bodies touch the surfaces meet perfectly with no voids or overlaps.

The first body is the Continental Crust. All onshore terrain is represented by the crust, and it also extends underneath ocean basins (see notes later). The top surface is the air-rock terrain boundary when onshore, and the water-rock bathymetry boundary when offshore. The total thickness of the crust when onshore is the sum of the terrain height (above, or in some cases below, sea level), plus a constant (35,000m), plus the isostatic root (or thinning) which is a function of the density contrast between the crust and mantle and the terrain height. When offshore, the thickness of the crust is given by a constant (35,000m), minus the bathymetric depth, minus the isostatic thinning which is a function of the density contrast between the crust, sea water, and the mantle and the bathymetric depth. The crust can include islands and other continental terrain lying offshore.

The diagram shown is from Wikipedia and illustrates the thickening and thinning of the crust due to terrain and bathymetry. The same Wikipedia page contains the formulas for the elevation vales b1 and b2 according to the Airy–Heiskanen model.

The second body is the Mantle. The top surface of the mantle is defined by the bottom surface of the crust, and the two surfaces mesh perfectly together. The bottom surface of the mantle is spherical at an arbitrary depth below sea level (200,000m).

The third body is the Ocean. The top surface of the ocean is sea level (0m RL). The bottom surface of the ocean is the water-rock bathymetry boundary. The ocean model is terminated at the shoreline. The ocean can contain islands and multiple enclosed basins.

The Continental Crust body extends into ocean basins where we would typically expect to find Oceanic Crust. However, I do not generate a body for oceanic crust. This can lead to an error as the continental crust density is considerably lower (2.83 g/cc) than oceanic crust (3.2 g/cc). To mitigate this issue, I use an accelerated crustal thinning formula that reduces the thickness of the continental crust to zero when offshore in deep ocean. The mantle body naturally replaces this missing crust and, with a density of 3.3 g/cc, is a more appropriate replacement for oceanic crust.

Input Variables and Raster Datasets

Density of continental crust: 2.83 g/cc
Density of mantle: 3.3 g/cc
Density of sea water: 1.025 g/cc
Thickness of continental crust: 35,000 m
Base of mantle: -200,000 m RL
Radius of model: 819,200 m
Complexity of model: 962,944 facets per body
Earth radius = 6371,008.7714 m
Deep ocean threshold: 4,000 m
Deep ocean thinning factor: 4X

Onshore terrain: MERIT_DEM

This is a near global 3 arc-second resolution raster that only includes onshore terrain and does not include Antarctica. I also used this raster as the Land-Water mask (to determine whether a negative terrain RL is onshore or offshore).

https://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_DEM/
https://robertsgeospatial.com.au/merit-dem-now-available-for-download-in-mrr-format/
https://robertsgeospatial.com.au/comparing-the-quality-cgiar-srtm-90m-dem-vs-merit-dem/

Offshore bathymetry: GEBCO 2023

This is a global 15 arc-second resolution raster that includes onshore terrain and offshore bathymetry.

https://www.gebco.net/

MERIT_DEM and GEBCO were combined into a single virtual raster which is then queried for terrain and bathymetry RL values. This virtual raster uses MERIT for onshore terrain, except below 59.5 degrees South (Antarctica) where it uses GEBCO. The virtual raster uses GEBCO for offshore bathymetry. In the near-shore environment where MERIT indicates the query location is offshore, but GEBCO returns a positive terrain value (due to its lower spatial resolution), the virtual raster returns -1 m RL.

The capture to the left shows this virtual raster illustrating high resolution onshore terrain, lower resolution offshore bathymetry, and an offshore discontinuity where an RL of -1 is returned to replace invalid bathymetry data.

Constructing the Earth Model

Results

The following four screen captures illustrate the results. All imagery is clipped to the NSW border. For a more regional context, see the map imagery at the top of the article.

The first map shows the onshore terrain.
The second map shows Free Air Gravity.
The third map shows Bouguer Anomaly (density unknown).
The fourth map shows the Corrected Gravity Anomaly. The Corrected Gravity Anomaly ranges from -58 to 71 milligals.

The Corrected Gravity Anomaly is computed for a subset (18.4%) of the total observation database. Stations were chosen based on their proximity to a 0.05-degree spaced grid. Even at this scale, you will see some difference in detail between the Bouguer Anomaly image and the Corrected Gravity image due to the difference in the number and density of observations gridded.

The model generated used a Crust density of 2.83 g/cc, Mantle density of 3.3 g/cc, and Ocean density of 1.025 g/cc. The model is based on a spherical approximation of the Earth. The correction, computed by calculating the vertical gravity component of the 3D faceted multi-body model at each observation location, combines the Bouguer, Terrain, and Isostatic corrections. The correction is subtracted from the Free Air Gravity at each observation to yield the Corrected Gravity Anomaly.

The isostatic correction has effectively removed the significant positive Bouguer anomaly that runs along the east coast. The negative Bouguer anomaly that tracks the Great Dividing Range has been turned into a positive anomaly in the main, perhaps attributable to basalts and mafic rocks throughout the Range. The Sydney Basin and Hunter Valley is negative or neutral reflecting the sandstone in this region.

The following four screen captures illustrate the vertical gravity component field computed from the multi-body model. All imagery is clipped to the NSW border. In all cases, I illustrate the raster data using a 16-color set of ranges with boundaries determined by the Jenks Natural Breaks algorithm.

The first map shows the field of the Crust model. The contribution ranges from 3809 to 4928 milligals.

The second map shows the field of the Ocean model. The contribution ranges from 0 to 9 milligals.

The third map shows the field of the Mantle model. The contribution ranges from 19620 to 20723 milligals.

The fourth map shows the field of the combined model. The gravity component ranges from 24422 to 24562 milligals.

Regional Traverses

I created four roughly West–East regional traverses. Along each traverse, I gathered nearby gravity survey locations on a 0.05-degree distance step. Thereafter, the station locations define the traverse. From south to north, the traverses are named Snowy Mountains, Blue Mountains, Central Tablelands, and Northern Rivers.

Along each traverse I computed the isostatic gravity correction for a crustal density of 2.67 and 2.83 g/cc and for a standard crustal thickness of 25, 30, 35, 40, 45, 50, 55, and 60 kilometres. The total thickness of the Crust at each gravity station is given by the standard thickness + terrain RL + Airy isostasy adjustment. A map of the traverses on top of terrain, clipped to the boundary of New South Wales, is shown below.

The four sections are illustrated below. The Gravity Anomaly has been computed for a crustal density of 2.83 g/cc. The brown curve is for a standard crustal thickness of 25 km and the green curve is for a standard crustal thickness of 55 km. The statistics of the computed gravity correction yields some insights.

Firstly, the standard deviation of the gravity is smaller when you use a crustal density of 2.83 vs 2.67 g/cc. This does not indicate what the most appropriate crustal density is, but I conclude it is closer to 2.83 than 2.67 g/cc.

Secondly, the standard deviation of the gravity decreases with increasing standard crustal thickness. I think this casts doubt on the usefulness of the Airy isostatic crustal variation away from the continental edge. On the continental edge, where it models a transition from continental crust to oceanic crust, it corrects a major residual Bouguer anomaly. However, in the continental interior, it seems to do little but introduce noise.

I have considered using a crustal thickness model based on geophysical data, AusMoho, to either augment or replace the Airy adjustments. The AusMoho model is shown as a contoured image overleaf with crustal thickness indicated in kilometres and contour interval of 0.5 kilometres. I think it is interesting that the AusMoho model shows little correlation with surface terrain (for example the thickening in the eastern Murray-Darling basin, the thickening in the Central West and New England North West, and the thinning along the North Coast).

Further Recommendations

This method is computationally intensive and improvements to the methodology ought to focus on reducing computation time. Building the model at each location is where efficiencies can be found.

It may be better to generate a single “super-model” for each UTM zone over which gravity observations are located. Although this model would be at least 10X larger and the time taken to compute the gravitational field of the model would take longer, the computation time overall would be reduced. In addition, any difference between the terrain raster and the gravity station RL could be incorporated into the model as a terrain perturbation.

We generate the model in projected UTM coordinates, but all the global raster coordinate systems are geographic. Consequently, all UTM coordinates are converted into geographic coordinates at each vertex of the model. It may be far more efficient to use virtual rasters to convert the rasters to the required projected coordinate system, eliminating the requirement to convert the vertex coordinates.

Using a crustal thickness model based on multiple geophysical observations (for example, AusMoho for Australia) may be a superior approach to using the Airy–Heiskanen approach.

This processing approach is useful for regional gravity surveys – like those undertaken by government research organisations – where isostatic effects can dominate the Bouguer Anomaly at large scales.