# gravfft

Spectral calculations of gravity, isostasy, admittance, and coherence for grids

## Synopsis

**gmt gravfft** *ingrid* [ *ingrid2* ]
**-G***outgrid*
[ **-C***n/wavelength/mean_depth*/**t**|**b**|**w** ]
[ **-D***density*|*rhogrid* ]
[ **-E***n_terms* ]
[ **-F**[**f**[**+s**|**z**]|**b**|**g**|**v**|**n**|**e**] ]
[ **-I****w**|**b**|**c**|**t**|**k** ]
[ **-N***params* ]
[ **-Q** ]
[ **-S** ]
[ **-T***te/rl/rm/rw*[*/ri*][**+m**] ]
[ **-V**[*level*] ]
[ **-W***wd*[**k**] ]
[ **-Z***zm*[*zl*] ]
[ **-f**flags ]
[ **--PAR**=*value* ]

**Note**: No space is allowed between the option flag and the associated arguments.

## Description

**gravfft** can be used into three main modes. Mode 1: Simply
compute the geopotential due to the surface given in the topo.grd file.
Requires a density contrast (**-D**) and possibly a different observation
level (**-W**). It will take the 2-D
forward FFT of the grid and
use the full Parker’s method up to the chosen terms. Mode 2: Compute the
geopotential response due to flexure of the topography file. It will take the 2-D
forward FFT of the grid and use the full Parker’s method applied to the chosen
isostatic model. The available
models are the “loading from top”, or elastic plate model, and the
“loading from below” which accounts for the plate’s response to a
sub-surface load (appropriate for hot spot modeling - if you believe
them). In both cases, the model parameters are set with **-T** and
**-Z** options. Mode 3: compute the admittance or coherence between
two grids. The output is the average in the radial direction.
Optionally, the model admittance may also be calculated.
Given the number of choices this program offers, is difficult to state
what are options and what are required arguments. It depends on what you
are doing; see the examples for further guidance.

## Required Arguments

*ingrid*[=*ID*|?*varname*][**+b***band*][**+d***divisor*][**+n***invalid*][**+o***offset*][**+s***scale*]

2-D binary grid file to be operated on. For cross-spectral operations, also give the second grid file

ingrid2. Optionally, append =IDfor reading a specific file format [Default is =nf] or ?varnamefor a specific netCDF variable [Default is the first 2-D grid found by GMT] (See full description). The following modifiers are supported:

+b- Select aband(for images only) [Default is 0].

+d- Divide data values by the givendivisor[Default is 1].

+n- Replace data values matchinginvalidwith NaN.

+o- Offset data values by the givenoffset[Default is 0].

+s- Scale data values by the givenscale[Default is 1].

Note: Any offset is added after any scaling.

**-G***outgrid*[=*ID*][**+d***divisor*][**+n***invalid*][**+o***offset*|**a**][**+s***scale*|**a**][:*driver*[*dataType*][**+c***options*]]

Specify the name of the output grid file or the 1-D spectrum table (see

-E) Optionally, append =IDfor writing a specific file format (See full description). The following modifiers are supported:

+d- Divide data values by givendivisor[Default is 1].

+n- Replace data values matchinginvalidwith a NaN.

+o- Offset data values by the givenoffset, or appendafor automatic range offset to preserve precision for integer grids [Default is 0].

+s- Scale data values by the givenscale, or appendafor automatic scaling to preserve precision for integer grids [Default is 1].

Note: Any offset is added before any scaling.+saalso sets+oa(unless overridden). To write specific formats via GDAL, use =gdand supplydriver(and optionallydataType) and/or one or more concatenated GDAL-cooptions using+c. See the “Writing grids and images” cookbook section for more details.

## Optional Arguments

**-C***n/wavelength/mean_depth*/**t**|**b**|**w**Compute only the theoretical admittance curves of the selected model and exit.

*n*and*wavelength*are used to compute (n * wavelength) the total profile length in meters.*mean_depth*is the mean water depth. Append data flags (one or two) of**t**|**b**|**w**in any order. Here,*t*= use “from top” model,*b*= use “from below” model. Optionally specify*w*to write wavelength instead of frequency.

**-D***density*|*rhogrid*Sets density contrast across surface. Used, for example, to compute the gravity attraction of the water layer that can later be combined with the free-air anomaly to get the Bouguer anomaly. In this case do not use

**-T**. It also implicitly sets**-N+h**. Alternatively, specify a co-registered grid with density contrasts if a variable density contrast is required.**Note**: Any NaNs found in the density grid will be replaced with the minimum density found.

**-E***n_terms*Number of terms used in Parker expansion (limit is 10, otherwise terms depending on n will blow out the program) [Default = 3].

**-F**[**f**[**+s**|**z**]|**b**|**g**|**v**|**n**|**e**]Specify desired geopotential field: compute geoid rather than gravity

**f**= Free-air anomalies (mGal) [Default]. Append**+s**to add in the slab implied when removing the mean value from the topography. This requires zero topography to mean no mass anomaly. Alternatively, to force the far-field to be exactly zero (i.e., the corner nodes of the grid), select**+z**instead.**b**= Bouguer gravity anomalies (mGal).**g**= Geoid anomalies (m).**v**= Vertical Gravity Gradient (VGG; 1 Eotvos = 0.1 mGal/km).**e**= East deflections of the vertical (micro-radian).**n**= North deflections of the vertical (micro-radian).

**-I****w**|**b**|**c**|**t**|**k**Use

*ingrid2*and*ingrid1*(a grid with topography/bathymetry) to estimate admittance|coherence and write it to standard output (**-G**ignored if set). This grid should contain gravity or geoid for the same region of*ingrid1*. Default computes admittance. Output contains 3 or 4 columns. Frequency (wavelength), admittance (coherence) one sigma error bar and, optionally, a theoretical admittance. Append dataflags (one to three) from**w**|**b**|**c**|**t**.**w**writes wavelength instead of wavenumber,**k**selects km for wavelength unit [m],**c**computes coherence instead of admittance,**b**writes a fourth column with “loading from below” theoretical admittance, and**t**writes a fourth column with “elastic plate” theoretical admittance.

**-N**[**a**|**f**|**m**|**r**|**s**|*nx/ny*][**+a**|**d**|**h**|**l**][**+e**|**n**|**m**][**+t***width*][**+v**][**+w**[*suffix*]][**+z**[**p**]]Choose or inquire about suitable grid dimensions for FFT and set optional parameters. Control the FFT dimension via these directives:

**a**: Let the FFT select dimensions yielding the most accurate result.**f**: Force the FFT to use the actual dimensions of the data.**m**: Let the FFT select dimensions using the least work memory.**r**: Let the FFT select dimensions yielding the most rapid calculation.**s**: Just present a list of optional dimensions, then exit.

Without a directive we expect

**-N***nx/ny*which will do FFT on array size*nx/ny*(must be >= grid file size). Default chooses dimensions >= data which optimize speed and accuracy of FFT. If FFT dimensions > grid file dimensions, data are extended and tapered to zero.Control detrending of data by appending a modifier for removing a linear trend. Consult module documentation for the default action:

**+d**: Detrend data, i.e. remove best-fitting linear trend.**+a**: Only remove the mean value.**+h**: Only remove the mid value, i.e. 0.5 * (max + min).**+l**: Leave data alone.

Control extension and tapering of data by appending a modifier to control how the extension and tapering are to be performed:

**+e**: Extend the grid by imposing edge-point symmetry [Default].**+m**: Extends the grid by imposing edge mirror symmetry.**+n**: Turns off data extension.

Tapering is performed from the data edge to the FFT grid edge [100%]. Change this percentage via modifier

**+t***width*. When**+n**is in effect, the tapering is applied instead to the data margins as no extension is available [0%].Control messages being reported:

**+v**: Report suitable dimensions during processing.

Control writing of temporary results: For detailed investigation you can write the intermediate grid being passed to the forward FFT; this is likely to have been detrended, extended by point-symmetry along all edges, and tapered. Use these modifiers to ave such grids:

**+w**: Set the*suffix*from which output file name(s) will be created (i.e.,*ingrid_prefix.ext*) [Default is “tapered”], where*ext*is your file extension**+z**: Save the complex grid produced by the forward FFT. By default we write the real and imaginary components to*ingrid*_real.*ext*and*ingrid*_imag.*ext*. Append**p**to instead use the polar form of magnitude and phase to files*ingrid*_mag.*ext*and*ingrid*_phase.*ext*.

**-Q**Writes out a grid with the flexural topography (with z positive up) whose average depth was set by

**-Z***zm*and model parameters by**-T**(and output by**-G**). That is the “gravimetric Moho”.**-Q**implicitly sets**-N+h**.

**-S**Computes predicted gravity or geoid grid due to a subplate load produced by the current bathymetry and the theoretical model. The necessary parameters are set within

**-T**and**-Z**options. The number of powers in Parker expansion is restricted to 1. See an example further down.

**-T***te/rl/rm/rw*[*/ri*][**+m**]Compute the isostatic compensation from the topography load (input grid file) on an elastic plate of thickness

*te*. Also append densities for load, mantle, water and infill in SI units. If*ri*is not provided it defaults to*rl*. Give average mantle depth via**-Z**. If the elastic thickness is > 1e10 it will be interpreted as the flexural rigidity (by default it is computed from*te*and Young modulus). Optionally, append*+m*to write a grid with the Moho’s geopotential effect (see**-F**) from model selected by**-T**. If*te*= 0 then the Airy response is returned.**-T+m**implicitly sets**-N+h**.

**-W***wd*[**k**]Set water depth (or observation height) relative to topography in meters [0]. Append

**k**to indicate km.

**-Z***zm*[/*zl*]Moho [and swell] average compensation depths (in meters positive down – the depth). For the “load from top” model you only have to provide

*zm*, but for the “loading from below” don’t forget to supply*zl*.

**-V**[*level*]Select verbosity level [

**w**]. (See full description) (See cookbook information).

**-f**flagsGeographic grids (dimensions of longitude, latitude) will be converted to meters via a “Flat Earth” approximation using the current ellipsoid parameters.

**-^**or just**-**Print a short message about the syntax of the command, then exit (

**Note**: on Windows just use**-**).**-+**or just**+**Print an extensive usage (help) message, including the explanation of any module-specific option (but not the GMT common options), then exit.

**-?**or no argumentsPrint a complete usage (help) message, including the explanation of all options, then exit.

**--PAR**=*value*Temporarily override a GMT default setting; repeatable. See gmt.conf for parameters.

## Grid Distance Units

If the grid does not have meter as the horizontal unit, append **+u***unit*
to the input file name to convert from the specified unit to meter. If your
grid is geographic, convert distances to meters by supplying **-f**flags instead.

## Considerations

netCDF COARDS grids will automatically be recognized as geographic. For
other grids geographical grids were you want to convert degrees into
meters, select **-f**flags. If the data are close to either pole, you should
consider projecting the grid file onto a rectangular coordinate system
using grdproject.

## Handling of Grids with NaNs

Since we cannot take FFTs of 2-D grids that contain NaNs, we perform simple substitutions.
If any of the input grids contain NaNs they will be replaced with zeros. In contrast, if **-D**
passes a grid with density contrasts then we replace any NaNs with the minimum density in the grid.

## Data Detrending

The default detrending mode follows Parker [1972] and removes the mid-value (**+h**).
Consult and use **-N** to select other modes.

## Plate Flexure

The FFT solution to elastic plate flexure requires the infill density to equal
the load density. This is typically only true directly beneath the load; beyond the load
the infill tends to be lower-density sediments or even water (or air). Wessel [2001]
proposed an approximation that allows for the specification of an infill density
different from the load density while still allowing for an FFT solution. Basically,
the plate flexure is solved for using the infill density as the effective load density but
the amplitudes are adjusted by a factor *A* = sqrt ((rm - ri)/(rm - rl)), which is
the theoretical difference in amplitude due to a point load using the two different
load densities. The approximation is very good but breaks down for large
loads on weak plates, a fairy uncommon situation.

## Examples

To compute the effect of the water layer above the bat.grd bathymetry using 2700 and 1035 for the densities of crust and water and writing the result on water_g.grd (computing up to the fourth power of bathymetry in Parker expansion):

```
gmt gravfft bat.grd -D1665 -Gwater_g.grd -E4
```

Now subtract it from your free-air anomaly faa.grd and you will get the Bouguer anomaly. You may wonder why we are subtracting and not adding. After all the Bouguer anomaly pretends to correct the mass deficiency presented by the water layer, so we should add because water is less dense than the rocks below. The answer relies on the way gravity effects are computed by the Parker’s method and practical aspects of using the FFT.

```
gmt grdmath faa.grd water_g.grd SUB = bouguer.grd
```

Want an MBA anomaly? Well compute the crust mantle contribution and add
it to the sea-bottom anomaly. Assuming a 6 km thick crust of density
2700 and a mantle with 3300 density we could repeat the command used to
compute the water layer anomaly, using 600 (3300 - 2700) as the density
contrast. But we now have a problem because we need to know the mean
Moho depth. That is when the scale/offset that can be appended to the grid’s name
comes in hand. Notice that we didn’t need to do that before because mean water
depth was computed directly from data (notice also the negative sign of the
offset due to the fact that *z* is positive up):

```
gmt gravfft bat.grd=+o-6000 -D600 -Gmoho_g.grd
```

Now, subtract it from the Bouguer to obtain the MBA anomaly. That is:

```
gmt grdmath bouguer.grd moho_g.grd SUB = mba.grd
```

To compute the Moho gravity effect of an elastic plate bat.grd with Te = 7 km, density of 2700, over a mantle of density 3300, at an average depth of 9 km

```
gmt gravfft bat.grd -Gelastic.grd -T7000/2700/3300/1035+m -Z9000
```

If you add now the sea-bottom and Moho’s effects, you will get the full gravity response of your isostatic model. We will use here only the first term in Parker expansion.

```
gmt gravfft bat.grd -D1665 -Gwater_g.grd -E1
gmt gravfft bat.grd -Gelastic.grd -T7000/2700/3300/1035+m -Z9000 -E1
gmt grdmath water_g.grd elastic.grd ADD = model.grd
```

The same result can be obtained directly by the next command. However,
PAY ATTENTION to the following. I don’t yet know if it’s because of a
bug or due to some limitation, but the fact is that the following and
the previous commands only give the same result if **-E**1 is used.
For higher powers of bathymetry in Parker expansion,
only the above example seams to give the correct result.

```
gmt gravfft bat.grd -Gmodel.grd -T7000/2700/3300/1035 -Z9000 -E1
```

And what would be the geoid anomaly produced by a load at 50 km depth, below a region whose bathymetry is given by bat.grd, a Moho at 9 km depth and the same densities as before?

```
gmt gravfft topo.grd -Gswell_geoid.grd -T7000/2700/3300/1035 -Fg -Z9000/50000 -S -E1
```

To compute the admittance between the topo.grd bathymetry and faa.grd free-air anomaly grid using the elastic plate model of a crust of 6 km mean thickness with 10 km effective elastic thickness in a region of 3 km mean water depth:

```
gmt gravfft topo.grd faa.grd -It -T10000/2700/3300/1035 -Z9000
```

To compute the admittance between the topo.grd bathymetry and geoid.grd geoid grid with the “loading from below” (LFB) model with the same as above and sub-surface load at 40 km, but assuming now the grids are in geographic and we want wavelengths instead of frequency:

```
gmt gravfft topo.grd geoid.grd -Ibw -T10000/2700/3300/1035 -Z9000/40000 -fg
```

To compute the gravity theoretical admittance of a LFB along a 2000 km long profile using the same parameters as above

```
gmt gravfft -C400/5000/3000/b -T10000/2700/3300/1035 -Z9000/40000
```

## References

Luis, J.F. and M.C. Neves. 2006, The isostatic compensation of the Azores Plateau: a 3D admittance and coherence analysis. J. Geothermal Volc. Res. Volume 156, Issues 1-2, Pages 10–22, https://doi.org/10.1016/j.jvolgeores.2006.03.010

Parker, R. L., 1972, The rapid calculation of potential anomalies, Geophys. J., 31, 447–455.

Wessel. P., 2001, Global distribution of seamounts inferred from gridded Geosat/ERS-1 altimetry, J. Geophys. Res., 106(B9), 19,431–19,441, https://doi.org/10.1029/2000JB000083