Issue about the result of hm.tesseroid_gravity function

[Arin47k]

Hello sir,
I am a PhD student. I am using the hm.tesseroid_gravity function for my research. Theoretically. If the volume of the tesseroid is zero, the we should get 'zero' gravity value at any outsider point of the tesseroid. But in this code, I get some value for the tesseroid with zero volume (bottom=top). Why this is happening, I can't get it. If you kindly help me in this regard, I shall be very grateful to you.
Thank you.
Regards,
Arindam
The code is given below:

import boule as bl
import pygmt
import verde as vd

import harmonica as hm

# Use the WGS84 ellipsoid to obtain the mean Earth radius which we'll use to
# reference the tesseroid
ellipsoid = bl.WGS84
mean_radius = ellipsoid.mean_radius

# Define tesseroid with top surface at the mean Earth radius, a thickness of
# 10km and a linear density function
tesseroids = (
    [-70, -60, -40, -30, mean_radius, mean_radius],
    [-70, -60, -30, -20, mean_radius, mean_radius],
    [-60, -50, -40, -30, mean_radius, mean_radius],
    [-60, -50, -30, -20, mean_radius-12000, mean_radius],
)

# Define a linear density function. We should use the jit decorator so Numba
# can run the forward model efficiently.


density = 2670,2670,2670,2670



# Define computation points on a regular grid at 100km above the mean Earth
# radius
coordinates = vd.grid_coordinates(
    region=[-70, -50, -40, -20],
    shape=(2,2),
    extra_coords=100e3 + ellipsoid.mean_radius,
)

# Compute the radial component of the acceleration
gravity = hm.tesseroid_gravity(coordinates, tesseroids, density, field="g_z",radial_adaptive_discretization=False)
print(gravity)
grid = vd.make_xarray_grid(
    coordinates, gravity, data_names="gravity", extra_coords_names="extra"
)

# Plot the gravitational field
fig = pygmt.Figure()

title = "Downward component of gravitational acceleration"

with pygmt.config(FONT_TITLE="16p"):
    fig.grdimage(
        region=[-70, -50, -40, -20],
        projection="M-60/-30/10c",
        grid=grid.gravity,
        frame=["a", f"+t{title}"],
        cmap="viridis",
    )

fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"])

fig.coast(shorelines="1p,black")

fig.show()

The output values:
[[ 11.18975045 16.33402878]
[ 17.85224879 325.11204181]]

I hope I will hear from you soon.

[santisoler]

Hi @Arin47k . Thanks for opening this question. Sorry for the delayed reply, I was too busy in the last month and I took the last week off.

From what I see in your code you do have one tesseroid with non-zero volume (the last one in your list), therefore it's expected that you'll get a non-zero gravitational effect on every observation point.

If you remove that last tesseroid in your list and rerun your script, you'll actually get zero values in your result:

import boule as bl
import verde as vd
import harmonica as hm

ellipsoid = bl.WGS84
mean_radius = ellipsoid.mean_radius

tesseroids = (
    [-70, -60, -40, -30, mean_radius, mean_radius],
    [-70, -60, -30, -20, mean_radius, mean_radius],
    [-60, -50, -40, -30, mean_radius, mean_radius],
)

density = 2670,2670,2670

coordinates = vd.grid_coordinates(
    region=[-70, -50, -40, -20],
    shape=(2,2),
    extra_coords=100e3 + ellipsoid.mean_radius,
)

gravity = hm.tesseroid_gravity(coordinates, tesseroids, density, field="g_z",radial_adaptive_discretization=False)
print(gravity)

And the result would be zero for every observation point:

array([[0., 0.],
       [0., 0.]])

So, you're right that tesseroids with zero volume should not create any gravitational effect. And that's what the function is doing. The values you're getting are due to the single tesseroid with volume in your list.

Hope this answers your question.

PS: I edited your question so your code looks nicer 🙂