Equivalent Sources for Gridding & Filtering Free-air Gravity Disturbance Data

[lawrenceabird]

Hi there:

I am hoping to make use of Harmonica, and in particular the equivalent sources method to grid some free-air gravity disturbance data that I have. I was hoping to get some thoughts/advice on my current plan of attack for my use case:

1. I have free-air gravity disturbance data (referenced to the WGS84 ellipsoid). Assuming that after the free-air correction, the elevation of these data is 0 m, is it reasonable to grid these data from their flight lines onto a regular grid using equivalent sources, and assigning an elevation of 0 m to the observations?
2. For my use case, I then need to transfer these data to a consistent elevation above the topography. Is it appropriate to upward-continue these gridded data using the extra_coords variable in the vd.grid_coordinates function?
3. I would like to try removing the long-wave component of the gravity signal. Is it appropriate to use a similar approach, but push the equivalent sources to a deeper depth (e.g., 20 km) to isolate the longer wavelength component of the signal (and then remove this from the other field)?

I haven't worked with geophysical data before, so definitely appreciate any advice here!

Thanks.

[santisoler]

Hi @lawrenceabird! Thanks for opening this question.

Gridding gravity data is one of the main purposes of equivalent sources, so I'm sure you'll be able to use them.

Before I answer your questions I wanted to clarify some concepts. The free-air gravity disturbance doesn't actually exists: given raw gravity measurements you can calculate the gravity anomaly, the free-air gravity anomaly or the gravity disturbance.
By definition the gravity disturbance $\delta g(\mathbf{p})$ at an observation point $\mathbf{p}$ is the observed gravity $g_\text{obs}(\mathbf{p})$ minus the normal gravity $\gamma(\mathbf{p})$ at the same observation point:

$$ \delta g(\mathbf{p}) = g_\text{obs}(\mathbf{p}) - \gamma(\mathbf{p}) $$

While the gravity anomaly $\Delta g(\mathbf{p})$ is the difference between the observed gravity and the normal gravity $\gamma(\mathbf{q})$ evaluated at a point $\mathbf{q}$ that stands on the surface of the ellipsoid and has the same longitudinal and latitudinal coordinates as $\mathbf{p}$. Point $\mathbf{q}$ has zero geometric height (zero height above the ellipsoid).

$$ \Delta g(\mathbf{p}) = g_\text{obs}(\mathbf{p}) - \gamma(\mathbf{q}) $$

The gravity anomaly is usually corrected with the free-air correction term(s). Considering the first order free-air correction it can be computed as:

$$ \Delta g_\text{fa}(\mathbf{p}) = g_\text{obs}(\mathbf{p}) - \gamma(\mathbf{q}) + \alpha h$$

where $\alpha$ is a constant factor (0.3086 if gravity is in mGal) and $h$ is the geometric height of the observation point $\mathbf{p}$ (its height above the ellipsoid).
The free-air anomaly can be rewritten as:

$$ \Delta g_\text{fa}(\mathbf{p}) = g_\text{obs}(\mathbf{p}) - (\gamma(\mathbf{q}) - \alpha h) $$

and here we can easily see that the terms inside the parentheses are just a first order approximation of the normal gravity $\gamma(\mathbf{p})$:

$$ \gamma(\mathbf{p}) = \gamma(\mathbf{q}) + \frac{\partial \gamma}{\partial h} \Big|_\mathbf{q} h + O(h^2) \simeq \gamma(\mathbf{q}) - \alpha h $$

But the thing is that we have analytical solutions for $\gamma$ on any point outside the ellipsoid. Boule provides reference ellipsoids with the normal_gravity() method that allows us to compute it.

So, when you say:

Assuming that after the free-air correction, the elevation of these data is 0 m,

note that this assumption is not correct: applying the free-air correction doesn't "take down" the observation point to the ellipsoid's surface. The elevation of the observation points after the free-air correction is unchanged.

But answering your question...

is it reasonable to grid these data from their flight lines onto a regular grid using equivalent sources, and assigning an elevation of 0 m to the observations?

Yes, you can grid your free-air gravity anomalies values from their flight lines onto a regular grid, but you need to respect the observation heights.

For my use case, I then need to transfer these data to a consistent elevation above the topography. Is it appropriate to upward-continue these gridded data using the extra_coords variable in the vd.grid_coordinates function?

Yes! This is the way to upward continue data with equivalent sources. Remember that upward continuation is perfectly fine, but downward continuation should be avoided or taken with a lot of care. So set a grid height that is equal to the highest observation point or above.

I would like to try removing the long-wave component of the gravity signal. Is it appropriate to use a similar approach, but push the equivalent sources to a deeper depth (e.g., 20 km) to isolate the longer wavelength component of the signal (and then remove this from the other field)?

Totally. This is a good way to separate the regional from the residual field ensuring that both are still harmonic fields and the variation of the field due to changes in elevation of the observation points is preserved. In this tutorial we show a small example on how we use deep equivalent sources for that. It's is also recommended to use a high damping parameter along with a large depth.

Hope this helps!

[lawrenceabird]

Hi @santisoler! Thanks so much for your reply.

I'm a little confused about the difference between the free-air gravity anomaly and the gravity disturbance. Would these two not be synonymous at the reference ellipsoid? If the free-air correction corrects for the elevation of the observation, then is the free-air gravity anomaly not equivalent to the gravity disturbance on the ellipsoid surface?

The dataset that I have is several years old - it calculates the normal gravity on the ellipsoid (and not at observation heights). A free-air correction was then applied : (0.308768 - 0.000440 * sin^2(Latitude) - 0.000000144 * ObservationElevation) * ObservationElevation, and the "gravity disturbance" then calculated and reported. So, in this situation, would it not be appropriate to use 0 m as the elevation when gridding the data using equivalent sources since the values are no longer representative of the gravity disturbance at the Observation Elevation, but on the ellipsoid (due to the free-air correction)? The data are not levelled in any way, and the Observation Elevations vary by ~1000 km, so I assume that I would not get a good result using equivalent sources in this case as the model would not be able to fit the observations at such varying elevation?

The other bits make sense -- thanks for confirming those!

[santisoler]

Hi @lawrenceabird!

I just want to make a minor comment on the equations I detailed before. They correspond to the modern definition of the gravity anomaly. @leouieda pointed out the differences if you consider the classical definition of the anomaly in https://github.com/orgs/fatiando/discussions/84#discussioncomment-4350725

I'm a little confused about the difference between the free-air gravity anomaly and the gravity disturbance. Would these two not be synonymous at the reference ellipsoid?

If your observation point is located on the ellipsoid (zero geometric height) then the values of the gravity anomaly and the gravity disturbance at that particular point are equal. But since most of the gravity observations are taken above the topography, at flight lines (both above the ellipsoid) and at sea level (not at zero geometric height, but at zero orthonormal height) this fact turns out to be a curious case that cannot be translated in real world applications. Moreover, the fact that these two fields have the same value at a particular point doesn't make them equivalent: for example the gradient of the two around a point located on the ellipsoid are not necessarily equal.

If the free-air correction corrects for the elevation of the observation, then is the free-air gravity anomaly not equivalent to the gravity disturbance on the ellipsoid surface?

The free-air anomaly and the gravity disturbance do take the same value on points of the ellipsoid, but again they are not equivalent fields. Since the free-air anomaly is a Taylor series approximation of the normal gravity they are not equivalent fields.

The dataset that I have is several years old - it calculates the normal gravity on the ellipsoid (and not at observation heights). A free-air correction was then applied : (0.308768 - 0.000440 * sin^2(Latitude) - 0.000000144 * ObservationElevation) * ObservationElevation, and the "gravity disturbance" then calculated and reported. So, in this situation, would it not be appropriate to use 0 m as the elevation when gridding the data using equivalent sources since the values are no longer representative of the gravity disturbance at the Observation Elevation, but on the ellipsoid (due to the free-air correction)?

In the case you are working with the data have been corrected using a second order approximation of the normal gravity. Nevertheless, this term is just approximating how the normal gravity varies with height and latitude, but it's not "lowering" the observation point to the ellipsoid. So the way I look at it is that you shouldn't assign a zero height to those points: they were obtained at a non-zero height and they have suffered the removal of a second order approximation of the normal gravity, so they weren't upward or downward continued in any way yet.

The data are not levelled in any way, and the Observation Elevations vary by ~1000 km, so I assume that I would not get a good result using equivalent sources in this case as the model would not be able to fit the observations at such varying elevation?

Are you sure that the heights vary in ~1000km? If these were obtained in flight lines I would love to see that plane going up to 1000 km 🛩️🧑‍🚀 😁. I assume that you meant 1000 m (1km). If that's so, there's no problem on applying equivalent sources. It might be an issue if you are trying to work with different types of data: ground observations (very low heights) and satellite observations (heights of ~200km).

[santisoler]

If you want to know more about definitions of gravity anomaly and disturbance, Barthelmes (2013) is a nice resource.

[lawrenceabird]

Hey @santisoler! Thanks again for explaining this further -- I think I'm starting to understand!

Can the free-air correction be interpreted as adjusting normal gravity that is calculated on the ellipsoid UP to the measurement elevation (rather than adjusting the measured gravity DOWN to the ellipsoid)?

And yes, by bad -- flight lines differ by ~1000 m, not km! Great to hear you think equivalent sources is still a reasonable approach.

One question on removing the long-wave regional component - do you have a sense of whether it is better calculate the residuals on the points, before gridding, or grid both components and then calculate the residuals? I guess if the equivalent sources model generates a good prediction it shouldn't matter too much, but can see benefits of both approaches.

Thanks!

[santisoler]

Hi @lawrenceabird! Glad to be helpful!

Can the free-air correction be interpreted as adjusting normal gravity that is calculated on the ellipsoid UP to the measurement elevation (rather than adjusting the measured gravity DOWN to the ellipsoid)?

This interpretation is not factually wrong, but I would argue that it's not accurate. The free-air correction is adding a first (or second) order approximation terms to the normal gravity on the ellipsoid term. So in a sense it is adjusting it up to the observation height. But it's a approximated adjustment. The approximation doesn't hold if you are working with satellite heights for example.

One question on removing the long-wave regional component - do you have a sense of whether it is better calculate the residuals on the points, before gridding, or grid both components and then calculate the residuals?

I always recommend to grid your data as the last step in your process, and it should be carried out only if there's an actual need for it. Gridding always introduces errors and biases, it shouldn't be considered as equivalent to the data or taken very carefully if otherwise. So I would recommend you to fit the deep sources on the scattered points, predict their effect (the regional component) on the same scattered points and compute the residual as the difference between the data on those scattered points and the regional field on the same points. Once you have the residual field on the scattered points you can grid them for better visualization and interpretation.

Are you planning to carry out more steps after the gridding? What's the goal of gridding the data for your particular use case?