fitting_gravity.md

Fitting gravity anomalies

IGMAS+ is designed for analyzing the time-independent subsurface density contributions which cause the gravity field; hence, anomalies to be analyzed should already be corrected for the effects related to instrumental drifts, Earth tides, the “normal gravity” reflecting the effects of flattening and the centrifugal force as well as the free air and/or mass corrections.

???+ note "Note on magnetic modelling"

For magnetic field modeling it is important to eliminate the International Geomagnetic Reference Field (IGRF), which describes the Earth's main magnetic field generated in the Earth's core. In the following, some important remarks are provided that need to be considered when modelling Bouguer and Free Air gravity anomalies. Magnetic field modelling is performed like modelling of the Free Air anomalies, since normally no magnetic effects of the topographic masses are subtracted from the field measurements.

---

Complete Bouguer anomaly (CBA): the "simple" situation

1. The gravity effects of the topographical masses have been eliminated from the measured values by the mass correction. This means, the terrain correction is included, accounting for the deviation of local topographic features by a flat or a spherical slab. Therefore, the background of the topo-masses is drawn transparent and bright in Figure. Stations of the terrestrial measurements are indicated by the green triangles; satellite gravity is measured at the orbit level (blue stippled line in Figure).

2. The density model extends to the geoid/ellipsoid surface. The model stations lie on the topography, at the orbit height or at any other user defined height. It must be ensured that model stations are identical with the heights and positions of gravity stations in the field/airplane or satellite orbit. Note: The density model is bordered by a constant model surface (geoid) as well as by a constant model bottom surface if model is built by vertical sections (refer to chapter "Model bodies", approach (a)). In case where horizons are used for model building (chapter "Model bodies". approach (b)), we face a special situation which is described later: a body with the density “0” is added automatically by the software.

3. It is irrelevant whether the modelling of satellite gravity is done at the orbit height, or at an arbitrarily chosen level (black with green triangles in Figure) above the highest elevation. It should be noted, however, that in the case of downward continuation of the satellite gravity field (from the satellite height to the arbitrarily chosen level), the so-called omission error (geodetic term) increases the further the chosen level is moved downwards. Reason: the gravity at the orbit height does not contain any small gravity wavelengths anymore, so that errors/inaccuracies/etc. are increased when the field continues downwards (through the massless space) and overlay the increased measurement signal.

4. The reference density of a model can be set arbitrarily by the user if the model structure does not cause a significant boundary effect.

5. It becomes difficult if not all mass effects between the surface and the reference surface (geoid/ellipsoid) could be eliminated because the rock densities deviate from the correction density (usually 2670 km/m^3). This information cannot be derived directly from gravity field modelling but must be extracted from independent information (e.g., geological maps and/or rock density determinations). In this case these deviating volumes must be later remodelled with a differential density:

\Delta\rho = \rho_{\mathrm{rock~mass}}-\rho_{\mathrm{mass~correction}}

---

Free Air anomaly (FA): the "simple" situation

In contrast to CBA, in the "simple" case of Free Air anomaly no mass correction is performed, so all masses in the model are preserved - including the topographic masses. The surface of the density model is now topography-dependent.

The terrestrial stations remain in the positions and at the heights as shown in Figure. If original satellite gravity is used for modelling purposes, those remain at the orbit height. If a grid with satellite gravity values (e.g. from the ICGEM data bases) is used, it can be extended down to any user-defined level (black line in Figure). Again, be careful: the omission error must be considered.

???+ warning "Attention:" One can "design" the surface of the density model using the heights of the stations and subtracting 13 cm from each station height:

```math
(\mathrm{heights~of~gravity~stations}) – 13~\mathrm{cm} = (\mathrm{height~of~model~top~surface})
```

It is a TIME CONSUMING PROCEDURE to check for all model sections whether the surface of topographic masses is 13 cm below the stations. The IGMAS+ developers are currently working on a much simpler method to use 3D topographies in the modelling. The result will appear in one of the next releases.

The value of 13 cm comes from the height (above the "ground") of the gravity meter measuring system of LaCoste-Gravimeter. If gravimeters of other companies are used, the measuring system height must be modified. Otherwise, one can assume an "overall" model station height of 1 m above the model surface; the error in larger crust/lithosphere models is to be neglected.

???+ tip "Hint:"

If the user decides to follow the “layer” approach ((b) in [chapter "Model setup"](model_setup.md#model-bodies)) and load layers/horizons for building up an initial 3D density model, **IGMAS+** automatically closes the model body upwards with a constant surface (see **Top** in the Body Manager). This “Top” body has zero density to mimic the air masses; for more detailed information please refer to [this chapter](#the-top-body-in-the-density-model).

There are two possibilities to define the model station heights (see Figure):

• the user chooses a constant height directly representing the satellite orbit (grey and pale blue body together) or
• the user continues the satellite data to a height level which is individually chosen by the modeller with reference to the heights proposed in the table of chapter "Model stations".

The "Top" body can be eliminated manually in the IGMAS+ input file. However, we suggest staying with “0”-density body; it does not cause any effect on the modelled anomaly and will disappear in one of the next IGMAS+ releases.

... and still, if this body disturbs someone, one can just eliminate it by hand. One can delete its polygons on every section:

• go to the first section
• right click on "Top"
• then select "Remove"
• repeat this procedure on every section.

If the "0-Body" doesn't have a single polygon anymore, you can also delete it from the Body Manager, because then "Remove Body" option is no longer greyed out.

Reference density: in the case of modelling a Free Air anomaly, the reference density must be set to zero (i.e. equal to the density of the topmost model body representing air)! Hence, one should check if the lateral model extensions are appropriately chosen to minimize edge effects.

---

Free Air anomaly (FA): the "difficult" situation

We start from the same situation as we have already studied for the Bouguer and Free Air anomalies: stations are at the orbit level or on the terrestrial surface (see Figures for CBA (simple) and FA (simple)). BUT: the terrain surface is now taken from a digital elevation model (DEM) or an elevation grid available on the web.

The illustration shows that the heights of the measured stations do not always correspond to the heights (often averaged) of the terrain model. Thus, it can happen that stations are located within the masses formed by the DEM. If this grid is loaded as "layer/horizon" (compare previous example), then some of the stations are located inside the mass. This leads to errors.

A similar situation must be considered if the satellite field data were continued to a level below the maximum DEM topography height. To be absolutely sure that this will not happen, choose a level that is "guaranteed" to be above the highest DEM elevation. For support refer again to chapter "Model stations".

A remarkably similar problem can also occur if the model surface of water masses in an offshore modelling scenario is not exactly 0 meters and the modeling stations are not 13 cm above it.

???+ note "Hint" When using grids (a) for satellite gravity (e.g. EIGEN-6C4, GOCE) for comparison with model gravity and (b) for topography (e.g. MERIT, EMODnet, etc.), always make sure that the model stations are positioned exactly on the nodes of the elevation grid. It must be ensured that both grids (gravity and heights) have the same grid structure/geometry (same grid width, same projection etc.). The grid with the satellite gravity can be downloaded for any user-defined level (see Figures for CBA (simple), FA (simple) and FA (diffucult)) above the topography.

---

The "Top" body in the density model

Here we will explain the intentions of the "NULL body", which already played a role in previous chapters. This is only valid if layers/horizons were loaded. For this purpose, let us look at Figure. Originally, the IGMAS+ function was used to read in individual horizons. Figure shows a stack of horizons that were taken as output from other survey results, e.g., from an offshore seismic campaign:

The software sorts the stack of layers first with ⑦ (the lowest layer) to ① (the highest layer). Because IGMAS+ only processes closed bodies, it “closes” bodies between layers ⑦--⑥, ⑥--⑤, ⑤--④, etc. automatically with the consequence that the user has to input a layer (in this Figure the top of the “Top layer” is at the sea level) with constant user-defined height (Z-top = 0 m) above the layer ①. Likewise, a constant model bottom layer has to be defined below the layer ⑦ to close the model’s bottom. Both density-values are user-defined. The "Top" body automatically generated by IGMAS+ gets the density of the water layer: 1000 kg/m^3, the "Bottom" body -- an appropriate user-defined density \rho_z. In the case of a model that contains both onshore and offshore areas, this body, of course, must be defined accordingly.

The automatic completion of an offshore model with a "Top" layer can also be used to complete an onshore model "upwards" if we model a Free-Air anomaly. In this case, the "Top" layer has a density of 0 kg/m^3. Its upper limit is to be defined by the user and corresponds to the height at which the small green triangles lie or the blue stippled line is drawn (see Figures for CBA (simple), FA (simple) and FA (difficult)). These heights are always obligatory to be above the highest elevation of the topography. In this case of use of "interfaces/horizons" for the input of the model geometry, the "Top" layer in the model is reinterpreted as “topography layer”. IGMAS+ will then automatically form a body whose upper boundary ("Top" layer in Figure) is defined by the user; this body is given the density \rho = 0 kg/m^3 and therefore has no gravity effect on the stations lying on the topography (big green triangles in Figures for CBA (simple), FA (simple) and FA (difficult)).

---

Modifications of the density model

The calculated gravity field of a density model can be analyzed based on the respective residual gravity:

When changing the density model, either by changing the density or the geometry of model bodies, IGMAS+ automatically and instantaneously adjusts the calculated and residual gravity anomaly, which is displayed by a 3D or a 2D viewer (possibly together with additional constraining data). The 2D viewer always displays one of the working planes, which makes them the actual scenes of interactive geometrical modifications as implemented through:

1. moving, deleting or adding of vertices, or
2. dividing bodies through additional intersections.

Through its interactive mode, IGMAS+ is primarily designed for analyzing and adjusting the 3D density model by visual inspection of gravity anomalies. This clearly implies some level of subjectivity in the model evaluation but a major advantage in compression to automatic algorithms is that it gives the user more control and especially the ability to learn how different features influences the result. For this reason, in addition to changing the density model through a complicated try-and-error procedure, IGMAS+ also allows to invert for the density (of one or more bodies) by minimizing the residual of a model.

---

Remarks on the use of ICGEM gravity datasets

ICGEM stands for “International Centre for Global Earth Models” (Ince et al., 2019). For 15 years ICGEM is one of the five worldwide services coordinated by the International Gravity Field Service (IGFS) of the International Association of Geodesy (IAG). Static and temporal global gravity field models of the Earth are provided in a standardized format with a possibility to assign a DOI number and interactive calculation and visualization services of gravity field functionals are available.

For more information refer also to the instructive ICGEM poster presented at EGU-2019 and to the ICGEM documentation “Definition of functionals of the geopotential and their calculation from spherical harmonic models” (Barthelmes, 2013; see also Ince et al., 2019).

The online availability of global models opens many research possibilities worldwide. However, to use the data provided from global model grids correctly one should pay an extra attention for what dataset actually represents. In the following section we provide some hints regarding the use of ICGEM models for geophysical modelling.

The ICGEM documentation (Barthelmes, 2013) in a sophisticated mathematical-physical form shows how the different functionals and models of the gravity field are calculated. For understanding it is of great advantage to have certain geodetic knowledge. In the following we will try to give some useful hints in a very simplified manner. The first note refers to all users, who use the ICGEM anomaly gravity_anomaly_bg.

???+ note

The [ICGEM documentation](http://icgem.gfz-potsdam.de/str-0902-revised.pdf) provides the following comment for the calculation of the "simple Bouguer anomaly":

"The (simple) Bouguer gravity anomaly (Functional selection $\Longrightarrow$ `gravity_anomaly_bg`) is defined by the classical gravity anomaly minus the attraction of the Bouguer plate. Here it will be calculated by the spherical approximation of the classical gravity anomaly minus $2 \pi G \rho H$ (eqs. 107 and 126 of [STR09/02](http://icgem.gfz-potsdam.de/str-0902-revised.pdf)). The topographic heights $H(\lambda,\phi)$ are calculated from the spherical harmonic model of topography ([ETOPO1](https://www.ngdc.noaa.gov/mgg/global/)) used up to the same maximum degree as the gravity field model:

- For $H \ge 0$ (rock)  $\rightarrow$ $\rho = 2670$ kg/m$^3$,
- For $H < 0$ (water) $\rightarrow$ $\rho = (2670--1025)$ kg/m$^3$ is used.

The density contrast between ice and rock has not been taken into account $\Longrightarrow$ the results for Greenland and Antarctica are not correct."

Please have in mind, that H is the height of topography above the geoid!

Fortunately, the differences between geoidal and ellipsoidal heights on Earth are small, however, present everywhere. From the geophysical/gravimetric viewpoint, not only the disregard of the gravity effects of ice masses is incorrect, but also the calculation itself. Firstly, a "flat bouguer plate" is considered, not a spherical Bouguer plate. Secondly, no topographic correction is conducted, which can lead to considerable errors in mountainous areas (Andes, Alps, Apennine, Himalaya). On the other hand, the transition areas between the continental margins and the oceans are also subject to errors, as Figure shows.

This illustration shows that in ICGEM's gravity_anomaly_bg not all masses were correctly removed. Remnants remain in the topography over land (density 2670 kg/m^3) and at the sea (density 1645 kg/m^3). The latter residue can be minimized by modelling the sea water layer by a body with a density of 1645~kg/m^3.

Another remark refers to the meaning of the term anomaly in geodesy and in geophysics (see, e.g., Li & Götze, 2001). The illustration below makes the difference graphically clear. A geodetic anomaly always implies, that the gravity anomaly \Delta g is calculated at the geoidal surface which is NOT the "real" topographic surface (see Figure, left). This implies a downward continuation of g_P into the topographic masses which is incorrect -- at least from a geophysical point of view.

In contrast, geophysicists use the term anomaly (here, more precisely, Free Air anomaly) to determine the gravity difference \Delta g directly at the observation point P. For this purpose, it is necessary to perform an upward continuation (height h in Figure) from the ellipsoid into the measuring point level P. This is a mathematically-physically stable procedure and therefore is allowed.

???+ note The geophysical Free Air anomaly is the same in meaning as a geodetic disturbance.

So, if the model gravity is to be compared with the Free Air anomaly, in the ICGEM calculation service first the gravity model is selected under the heading "Model selection" (EIGEN-6C4 is a good choice, also XGM2019), then in the heading "Functional selection" disturbance is selected, if no specific non-zero meter altitude is to be selected.

This is dangerous, because then, under certain circumstances, the stations of the Free Air anomaly/disturbance can lie within the model masses (see chapters on FA (simple), FA (difficult) and on the "Top" body). If the model contains topographic heights which are higher than 0 m, then select disturbance_sa in the column "Functional selection" and enter the height in the input mask on the right under the input of "grid step [°]". The height is freely selectable by the user and must be above the maximum of the model topography; this corresponds to the height with the small green triangles in Figures for CBA (simple), FA (simple) and FA (difficult).

To calculate complete Bouguer anomalies it is recommended to download the free-air/disturbance at station level and then perform Bouguer corrections using a DEM model for both on- and offshore masses.