talwani2d

Compute geopotential anomalies over 2-D bodies by the method of Talwani

Synopsis

gmt talwani2d [ table ] [ -A ] [ -Ddensity ] ] [ -Ff|n[lat]|v ] [ -M[h][v] ] [ -Ntrackfile ] [ -Tmin/max/inc[+i|n]|file|list ] [ -Zlevel[/ymin/ymax] ] [ -V[level] ] [ -bibinary ] [ -dnodata[+ccol] ] [ -eregexp ] [ -hheaders ] [ -iflags ] [ -oflags ] [ -x[[-]n] ] [ --PAR=value ]

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

Description

talwani2d will read the multi-segment table from file (or standard input). This file contains cross-sections of one or more 2-D bodies, with one polygon per segment. The segment header must contain the parameter density, which states the density of this body (individual body densities may be overridden by a fixed constant density contrast given via an optional -D). We can compute anomalies on an equidistant lattice (by specifying a lattice with -T) or provide arbitrary output points specified in a file via -N. Choose between free-air anomalies, vertical gravity gradient anomalies, or geoid anomalies. Options are available to control axes units and direction. For theory, see references at the end.

Required Arguments

table: One or more ASCII files describing cross-sectional polygons of one or more bodies. Polygons will be automatically closed if not already closed, and repeated vertices will be eliminated. The segment header for each body will be examined for a density parameter in \(\mbox{kg/m}^3\) or \(\mbox{g/cm}^3\); see -D for overriding this value. If no table is given then we read standard input.

Optional Arguments

-A: The z-axis should be positive upwards [Default is down].

-Ddensity: Sets a fixed density contrast that overrides any per-body settings in the model file, in \(\mbox{kg/m}^3\) or \(\mbox{g/cm}^3\).

-Ff|n[lat]|v: Specify desired gravitational field component. Choose between f (free-air anomaly) [Default], n (geoid; optionally append average latitude for normal gravity reference value [45]) or v (vertical gravity gradient).

-M[h][v]: Sets distance units used. Append h to indicate horizontal distances are in km [m], and append z to indicate vertical distances are in km [m].

-Ntrackfile: Specifies locations where we wish to compute the predicted value. When this option is used you cannot use -T to set an equidistant lattice. The output data records are written to standard output (see -bo for binary output).

-Tmin/max/inc[+i|n]|file|list: Specify an equidistant output lattice. For details on array creation, see Generate 1-D Array.

-V[level]: Select verbosity level [w]. (See full description) (See technical reference).

-Zlevel[/ymin/ymax]: Set a constant observation level [0]. Optionally, and for gravity anomalies only (-Ff), append the finite extent limits of a 2.5-D body.

-birecord[+b|l] (more …): Select native binary format for primary table input. [Default is 2 input columns].

-d[i|o][+ccol]nodata (more …): Replace input columns that equal nodata with NaN and do the reverse on output.

-e[~]“pattern” | -e[~]/regexp/[i] (more …): Only accept data records that match the given pattern.

-h[i|o][n][+c][+d][+msegheader][+rremark][+ttitle] (more …): Skip or produce header record(s). Not used with binary data.

-icols[+l][+ddivisor][+sscale|d|k][+ooffset][,…][,t[word]] (more …): Select input columns and transformations (0 is first column, t is trailing text, append word to read one word only).

-ocols[+l][+ddivisor][+sscale|d|k][+ooffset][,…][,t[word]] (more …): Select output columns and transformations (0 is first column, t is trailing text, append word to write one word only).

-x[[-]n] (more …): Limit number of cores used in multi-threaded algorithms. Multi-threaded behavior is enabled by default. That covers the modules that implement the OpenMP API (required at compiling stage) and GThreads (Glib) for the grdfilter module.

-:[i|o] (more …): Swap 1st and 2nd column on input and/or output.

-^ 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 arguments: Print 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.

Units

For map distance unit, append unit d for arc degree, m for arc minute, and s for arc second, or e for meter [Default unless stated otherwise], f for foot, k for km, M for statute mile, n for nautical mile, and u for US survey foot. By default we compute such distances using a spherical approximation with great circles (-jg) using the authalic radius (see PROJ_MEAN_RADIUS). You can use -jf to perform “Flat Earth” calculations (quicker but less accurate) or -je to perform exact geodesic calculations (slower but more accurate; see PROJ_GEODESIC for method used).

Generate 1-D Array

We will demonstrate the use of options for creating 1-D arrays via math. Make an evenly spaced coordinate array from min to max in steps of inc, e.g.:

gmt math -o0 -T3.1/4.2/0.1 T =
1
2
3
4
5
6
7
...

Append +b if we should take \(\log_2\) of min and max, get their nearest integers, build an equidistant \(\log_2\)-array using inc integer increments in \(\log_2\), then undo the \(\log_2\) conversion. E.g., -T3/20/1+b will produce this sequence:

gmt math -o0 -T3/20/1+b T =
4
8
16

Append +l if we should take \(\log_{10}\) of min and max and build an array where inc can be 1 (every magnitude), 2, (1, 2, 5 times magnitude) or 3 (1-9 times magnitude). E.g., -T7/135/2+l will produce this sequence:

gmt math -o0 -T7/135/2+l T =
10
20
50
100

For output values less frequently than every magnitude, use a negative integer inc:

gmt math -o0 -T1e-4/1e4/-2+l T =
0.0001
0.01
1
100
10000

Append +i if inc is a fractional number and it is cleaner to give its reciprocal value instead. To set up times for a 24-frames per second animation lasting 1 minute, run:

gmt math -o0 -T0/60/24+i T =
0
0.0416666666667
0.0833333333333
0.125
0.166666666667
...

Append +n if inc is meant to indicate the number of equidistant coordinates instead. To have exactly 5 equidistant values from 3.44 and 7.82, run:

gmt math -o0 -T3.44/7.82/5+n T =
44
535
63
725
82

Alternatively, let inc be a file with output coordinates in the first column, or provide a comma-separated list of specific coordinates, such as the first 6 Fibonacci numbers:

gmt math -o0 -T0,1,1,2,3,5 T =
0
1
1
2
3
5

Notes: (1) If you need to pass the list nodes via a dataset file yet be understood as a list (i.e., no interpolation), then you must set the file header to contain the string “LIST”. (2) Should you need to ensure that the coordinates are unique and sorted (in case the file or list are not sorted or have duplicates) then supply the +u modifier.

If you only want a single value then you must append a comma to distinguish the list from the setting of an increment.

If the module allows you to set up an absolute time series, append a valid time unit from the list year, month, day, hour, minute, and second to the given increment; add +t to ensure time column (or use -f). Note: The internal time unit is still controlled independently by TIME_UNIT. The first 7 days of March 2020:

gmt math -o0 -T2020-03-01T/2020-03-07T/1d T =
2020-03-01T00:00:00
2020-03-02T00:00:00
2020-03-03T00:00:00
2020-03-04T00:00:00
2020-03-05T00:00:00
2020-03-06T00:00:00
2020-03-07T00:00:00

A few modules allow for +a which will paste the coordinate array to the output table.

Likewise, if the module allows you to set up a spatial distance series (with distances computed from the first two data columns), specify a new increment as inc with a geospatial distance unit from the list degree (arc), minute (arc), second (arc), meter, foot, kilometer, Miles (statute), nautical miles, or survey foot; see -j for calculation mode. To interpolate Cartesian distances instead, you must use the special unit c.

Finally, if you are only providing an increment and will obtain min and max from the data, then it is possible (max - min)/inc is not an integer, as required. If so, then inc will be adjusted to fit the range. Alternatively, append +e to keep inc exact and adjust max instead (keeping min fixed).

Examples

To compute the free-air anomalies on an equidistant profile over a 2-D body that has been contoured and saved to body2d.txt, using 1700 \(\mbox{kg/m}^3\) as a constant density contrast, with all distances in meters, try

gmt talwani2d -T-200/200/2 body2d.txt -D1700 -Ff > 2dgrav.txt

To obtain the vertical gravity gradient anomaly along the track given by the file crossing.txt for the same model, try

gmt talwani2d -Ncrossing.txt body2d.txt -D1700 -Fv > vgg_crossing.txt

The geoid anomaly for the same setup, evaluated at 60N, is given by

gmt talwani2d -Ncrossing.txt body2d.txt -D1700 -Fn60 > n_crossing.txt

Notes

The 2-D geoid anomaly is a logarithmic potential and thus has no natural reference level. We simply remove the most negative (if density contrast is positive) or positive (if density contrast is negative) computed value from all values, rendering the entire anomaly positive (or negative). You can use gmtmath to change the zero level to suit your needs.

References

Rasmussen, R., and L. B. Pedersen (1979), End corrections in potential field modeling, Geophys. Prospect., 27, 749-760.

Chapman, M. E., 1979, Techniques for interpretation of geoid anomalies, J. Geophys. Res., 84(B8), 3793-3801.

Kim, S.-S., and P. Wessel, 2016, New analytic solutions for modeling vertical gravity gradient anomalies, Geochem. Geophys. Geosyst., 17, https://doi.org/10.1002/2016GC006263.

Talwani, M., J. L. Worzel, and M. Landisman, 1959, Rapid gravity computations for two-dimensional bodies with application to the Mendocino submarine fracture zone, J. Geophys. Res., 64, 49-59.