Ctrl+K

harmonica.prism_gravity

harmonica.prism_gravity#

harmonica.prism_gravity(coordinates, prisms, density, field, parallel=True, dtype='float64', progressbar=False, disable_checks=False)[source]#

Gravitational fields of right-rectangular prisms in Cartesian coordinates

Compute the gravitational potential, gravitational acceleration and tensor components generated by a collection of prisms on a set of observation points.

Warning

The vertical direction points upwards, i.e. positive and negative values of upward represent points above and below the surface, respectively. But g_z field returns the downward component of the gravitational acceleration so that positive density contrasts produce positive anomalies. The same applies to the tensor components, i.e. the g_ez is the non-diagonal easting-downward tensor component.

Important

  • The gravitational potential is returned in \(\text{J}/\text{kg}\).

  • The gravity acceleration components are returned in mGal (\(\text{m}/\text{s}^2\)).

  • The tensor components are returned in Eotvos (\(\text{s}^{-2}\)).

Parameters:

  • coordinates (list of arrays) – List of arrays containing the easting, northing and upward coordinates of the computation points defined on a Cartesian coordinate system. All coordinates should be in meters.

  • prisms (list, 1d-array, or 2d-array) – List or array containing the coordinates of the prism(s) in the following order: west, east, south, north, bottom, top, in a Cartesian coordinate system. All coordinates should be in meters. Coordinates for more than one prism can be provided. In this case, prisms should be a list of lists or 2d-array (with one prism per row).

  • density (list or array) – List or array containing the density of each prism in kg/m^3.

  • field (str) –

    Gravitational field that wants to be computed. The available fields are:

    • Gravitational potential: potential

    • Eastward acceleration: g_e

    • Northward acceleration: g_n

    • Downward acceleration: g_z

    • Diagonal tensor components: g_ee, g_nn, g_zz

    • Non-diagonal tensor components: g_en, g_ez, g_nz

  • parallel (bool (optional)) – If True the computations will run in parallel using Numba built-in parallelization. If False, the forward model will run on a single core. Might be useful to disable parallelization if the forward model is run by an already parallelized workflow. Default to True.

  • dtype (data-type (optional)) – Data type assigned to the resulting gravitational field. Default to np.float64.

  • progressbar (bool (optional)) – If True, a progress bar of the computation will be printed to standard error (stderr). Requires numba_progress to be installed. Default to False.

  • disable_checks (bool (optional)) – Flag that controls whether to perform a sanity check on the model. Should be set to True only when it is certain that the input model is valid and it does not need to be checked. Default to False.

Returns:

result (array) – Gravitational field generated by the prisms on the computation points. Gravitational potential is returned in \(\text{J}/\text{kg}\), acceleration components in mGal, and tensor components in Eotvos.

Notes

This function makes use of choclo.prism forward modelling functions to compute each gravitational field.

The gravitational potential ("potential") and the acceleration components ("g_e", "g_n" and "g_z") are well defined on the entire domain. Tensor components aren’t defined on prism vertices. Diagonal tensor components aren’t defined on edges normal to the direction of the tensor (e.g. “g_ee” is not defined on edges parallel to northing and upward directions). Non-diagonal tensor components aren’t defined on edges normal to the remaining direction of the tensor (e.g. “g_en” is not defined on edges parallel to the upward direction). The function returns numpy.nan on every singular point.

The diagonal tensor components aren’t defined on observation points that belong to the faces normal to their direction (e.g. "g_zz" is not define on horizontal faces): two different limits exist when approaching from the inside and from the outside of the prism. This function returns the limit of these components while approaching from the outside.

References

Examples

Compute gravitational effect of a single a prism

>>> # Define prisms boundaries, it must be beneath the surface
>>> prism = [-34, 5, -18, 14, -345, -146]
>>> # Set prism density to 2670 kg/m³
>>> density = 2670
>>> # Define three computation points along the easting direction at 30m
>>> # above the surface
>>> easting = [-40, 0, 40]
>>> northing = [0, 0, 0]
>>> upward = [30, 30, 30]
>>> coordinates = (easting, northing, upward)
>>> # Compute the downward component of the gravitational acceleration that
>>> # the prism generates on the computation points
>>> gz = prism_gravity(coordinates, prism, density, field="g_z")
>>> print("({:.5f}, {:.5f}, {:.5f})".format(*gz))
(0.06552, 0.06629, 0.06174)

Define two prisms with positive and negative density contrasts

>>> prisms = [[-134, -5, -45, 45, -200, -50], [5, 134, -45, 45, -180, -30]]
>>> densities = [-300, 300]
>>> # Compute the g_z that the prisms generate on the computation points
>>> gz = prism_gravity(coordinates, prisms, densities, field="g_z")
>>> print("({:.5f}, {:.5f}, {:.5f})".format(*gz))
(-0.05380, 0.02908, 0.11237)

Examples using harmonica.prism_gravity#

Layer of prisms

Layer of prisms

Gravitational effect of topography

Gravitational effect of topography
Contents