Provided by: gmt-common_5.2.1+dfsg-3build1_all 
NAME
gravfft - Compute gravitational attraction of 3-D surfaces in the wavenumber (or frequency) domain
SYNOPSIS
gravfft ingrid [ ingrid2 ] outfile [ n/wavelength/mean_depth/tbw ] [ density|rhogrid ] [ n_terms ] [
[f[+]|g|v|n|e] ] [ w|b|c|t |k ] [ [f|q|s|nx/ny][+a|d|h |l][+e|n|m][+twidth][+w[suffix]][+z[p]] [ ] [
te/rl/rm/rw[+m] ] [ [level] ] [ wd] [ zm[zl] ] [ -fg ]
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. The horizontal dimensions of the grdfiles are
assumed to be in meters. Geographical grids may be used by specifying the -fg option that scales degrees
to meters. If you have grids with dimensions in km, you could change this to meters using grdedit or
scale the output with grdmath. 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 2-D binary grid file to be operated on. (See GRID FILE FORMATS below). For cross-spectral
operations, also give the second grid file ingrd2.
-Goutfile
Specify the name of the output grid file or the 1-D spectrum table (see -E). (See GRID FILE
FORMATS below).
OPTIONAL ARGUMENTS
-Cn/wavelength/mean_depth/tbw
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 dataflags (one or two) of tbw in any order. t = use "from top" model, b = use
"from below" model. Optionally specify w to write wavelength instead of frequency.
-Ddensity|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.
-En_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[+]|g|v|n|e]
Specify desired geopotential field: compute geoid rather than gravity
f = Free-air anomalies (mGal) [Default]. Append + to add in the slab implied when removing the
mean value from the topography. This requires zero topography to mean no mass anomaly.
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).
-Iw|b|c|t |k
Use ingrd2 and ingrd1 (a grid with topography/bathymetry) to estimate admittance|coherence and
write it to stdout (-G ignored if set). This grid should contain gravity or geoid for the same
region of ingrd1. 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[f|q|s|nx/ny][+a|[+d|h|l][+e|n|m][+twidth][+w[suffix]][+z[p]]
Choose or inquire about suitable grid dimensions for FFT and set optional parameters. Control the
FFT dimension:
-Nf will force the FFT to use the actual dimensions of the data.
-Nq will inQuire about more suitable dimensions, report those, then continue.
-Ns will present a list of optional dimensions, then exit.
-Nnx/ny 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: Append modifiers for removing a linear trend:
+d: Detrend data, i.e. remove best-fitting linear trend [Default].
+a: Only remove mean value.
+h: Only remove mid value, i.e. 0.5 * (max + min).
+l: Leave data alone.
Control extension and tapering of data: Use modifiers to control how the extension and tapering
are to be performed:
+e extends 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 +twidth. When +n is in effect, the tapering is applied instead to the data margins as no
extension is available [0%].
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. Append +w[suffix] from which output file name(s) will
be created (i.e., ingrid_prefix.ext) [tapered], where ext is your file extension. Finally, you may
save the complex grid produced by the forward FFT by appending +z. By default we write the real
and imaginary components to ingrid_real.ext and ingrid_imag.ext. Append p to save instead 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 was set by -Zzm
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.
-Tte/rl/rm/rw[+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, and water in SI units. 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
-Wwd Set water depth (or observation height) relative to topography [0]. Append k to indicate km.
-Zzm[zl]
Moho [and swell] average compensation depths. For the "load from top" model you only have to
provide zm, but for the "loading from below" don't forget zl.
-V[level] (more ...)
Select verbosity level [c].
-fg Geographic 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 exits (NOTE: on Windows use just -).
-+ or just +
Print an extensive usage (help) message, including the explanation of any module-specific option
(but not the GMT common options), then exits.
-? or no arguments
Print a complete usage (help) message, including the explanation of options, then exits.
--version
Print GMT version and exit.
--show-datadir
Print full path to GMT share directory and exit.
GRID FILE FORMATS
By default GMT writes out grid as single precision floats in a COARDS-complaint netCDF file format.
However, GMT is able to produce grid files in many other commonly used grid file formats and also
facilitates so called "packing" of grids, writing out floating point data as 1- or 2-byte integers. To
specify the precision, scale and offset, the user should add the suffix =id[/scale/offset[/nan]], where
id is a two-letter identifier of the grid type and precision, and scale and offset are optional scale
factor and offset to be applied to all grid values, and nan is the value used to indicate missing data.
In case the two characters id is not provided, as in =/scale than a id=nf is assumed. When reading
grids, the format is generally automatically recognized. If not, the same suffix can be added to input
grid file names. See grdconvert and Section grid-file-format of the GMT Technical Reference and Cookbook
for more information.
When reading a netCDF file that contains multiple grids, GMT will read, by default, the first
2-dimensional grid that can find in that file. To coax GMT into reading another multi-dimensional
variable in the grid file, append ?varname to the file name, where varname is the name of the variable.
Note that you may need to escape the special meaning of ? in your shell program by putting a backslash in
front of it, or by placing the filename and suffix between quotes or double quotes. The ?varname suffix
can also be used for output grids to specify a variable name different from the default: "z". See
grdconvert and Sections modifiers-for-CF and grid-file-format of the GMT Technical Reference and Cookbook
for more information, particularly on how to read splices of 3-, 4-, or 5-dimensional grids.
GRID DISTANCE UNITS
If the grid does not have meter as the horizontal unit, append +uunit to the input file name to convert
from the specified unit to meter. If your grid is geographic, convert distances to meters by supplying
-fg 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 -fg. If the data are close to either pole, you
should consider projecting the grid file onto a rectangular coordinate system using grdproject.
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 to 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=nf/1/-6000 -D600 -Gmoho_g.grd
Now, subtract it to the sea-bottom anomaly to obtain the MBA anomaly. That is:
gmt grdmath water_g.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 -E1 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 the 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,
http://dx.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,
http://dx.doi.org/10.1029/2000JB000083
SEE ALSO
gmt, grdfft, grdmath, grdproject
COPYRIGHT
2015, P. Wessel, W. H. F. Smith, R. Scharroo, J. Luis, and F. Wobbe