Terrain correction - prisms vs tesseroids layer

[elisekazmierczak-geo]

I’m currently producing a new Bouguer anomaly map for Belgium. To handle the terrain correction, I would like to compare the prism layer and tesseroid leyer methods using the same digital terrain model. All the gravimetric data I have were measured on the ground.

The issue I’m facing is that certain points are considered to be “inside” the tesseroids, which does not seem to be a problem when using prisms. What solution would you recommend to ensure a coherent comparison?

My current approach is to assign each point the height of the surface of the tesseroid in which it is located. However, this means I’m using different heights in the two methods. Wouldn’t it be better, then, to apply the same procedure to the prism corrections—by assigning each measurement point the height of the prism surface?

Although the heights would still differ from those used with tesseroids, I assume that the comparison between the two terrain correction methods would then be based solely on the digital terrain model, making it more consistent.

Thank you for your help!

Elise

[santisoler]

Hi Elise. Thanks for raising this question!

Some background first. The gravity fields of prisms are computed using an analytic solution that is valid inside and outside of them. The gravity field of tesseroids on the other hand are computed using numerical approximations of the integrals. Those integrals are not accurately tested on internal points, and that's why we advise against computing them on internal points.

I think a potential source of the issue you are experiencing is that the digital elevation model (DEM) is not 100% in agreement with the height measurement when the gravity survey was performed: apparently some observation points are below the topographic height provided by the DEM. I'd first double check that the DEM and the heights of the observation points are defined in the same way (e.g. ellipsoidal heights vs geodetic heights). If they both are ellipsoidal heights, and the differences between them are not significant, then your approach of correcting either the height of the observation point or the height of the DEM for that point sounds reasonable.

Applying the same technique to both the prism and the tesseroid models seems reasonable too. Basically you are correcting for differences between the two measured heights (the ones of the DEM and the ones of the observation points) and then building models to compute the terrain correction on every observation point.

Hope it helps! And please let us know if you have any other question.

[elisekazmierczak-geo]

Hi Santiago,

Thank you very much for taking the time to reply to me!
All these clarifications are extremely helpful for my project.

The elevation system I use is consistent for both the DEM and the observation points. The gravimetric data was pre-filtered when the elevation differences with the DEM were too large (I used an approach similar to that of Medved et al., 2021). However, some points are still located slightly below the DEM surface. Additionally, many points lie in valleys, and when I reduce the resolution of my DEM to discretize the topography into tesseroids, the interpolation places the surface above the actual terrain level. This resolution reduction is also dependent on the extent of my DEM, which brings me to my second question.

If I understood correctly, neither the prism nor tesseroid layer methods apply any special treatment to the area closest to the observation point and the prism-layer method doesn’t account for Earth curvature effects. I’ve tested several extension distances around the calculation point (and, consequently, different DEM extents), as the literature seems to vary quite a bit (e.g., 167 km, 22 km, etc.), and these choices can lead to differences on the order of several mGal.
Do you have a recommended extension distance specific to these methods?

Thanks again for your time. I’d be delighted if my research could contribute to the Fatiando project—perhaps by publishing my codes and data as examples.

Cheers from Brussels,
Elise

[craigmillernz]

Hi Elise
A couple of suggestions to try. To deal with the issue of the station height and DEM not matching (assuming same datums etc), then i usually extract the DEM height at the station location and use that as the station elevation for the terrain correction. This ensures you station isn't "inside" the topography or floating a long way above it. You can add a small offset if you like (0.1 m) to account for your meter being a little bit off the ground. They key though is to use a fine enough resolution DEM so that the dem cell value at the station location isn't a large average of the surrounding terrain. In NZ we typically use a lidar grid (-.5 to 2m resolution) for our "inner zones" out to 170 m and I replace the measured station elevation with the DEM value at the station location using those fine grids. Then use a DEM of larger cell size (8m in NZ) out to 2610 m, then another grid (200 m) out to 21.9 km and finally a 2km DEM out to 167 km. Only clamp the station elevation to the DEM using the finest grid, and use that elevation for all other grids. 167 km is equivalent to 1.5 degrees, 21.9 km is "Hammer zone M" and 2.61 km is "Hammer zone H", while 170m is "Hammer zone D" in the traditional set up. Of course these distances are arbitary, but the Hammer zones are useful for reference. Will you include bathymetry in your terrain corrections?

[elisekazmierczak-geo]

Hi Craig, Thank you so much for your comprehensive and accurate response, it is extremely helpful! For the moment, I am using a coarse DEM with a resolution of several tens of meters that does not take bathymetry into account (fortunately, Belgium only has 65 km of coastline!). The idea is to first carry out a low resolution study and test several topographic correction methods before moving on to finer resolution work, which will be more computationally intensive. If I understand your method correctly, you cut your finest DEM into a circle centered on the measurement point and then you ‘cut’ your other coarser DEMs into rings around the disc centered on the measurement point. You then calculate the topographic correction for each circle/ring you are using and then add them together? Have I understood correctly? Do you use the Fatiando library to do this, or do you do it from another program? I will keep all your good advice in mind for the rest of my work, thank you again! :) Cheers from Brussels!

[craigmillernz]

Hi Elise, if you just have a coarse grid to start with then you can cut the DEM to the full distance and just use one "ring". Otherwise yes, it's using nested rings. This is mainly for computation requirements to reduce the amount of memory needed, but also for practicality as the you don't need to model every topographic detail once you are further from the station. As the grids are read in as xarrays i just functions in xarray to cut the DEM to size around each data point using the grid and station coordinates and simple geometry to calculate the new grid extents around each station. You can do "square" or "circular" girds. Note that if you set the base of the prism to 0 m and the top to the DEM elevation then you will calculate the full "topographic" correction which means you don't need the Bouguer slab correction in your processing.

[elisekazmierczak-geo]

Hi Craig! Thank you very much for all these valuable details. They help me see things more clearly! I was already using rioxarray functions to cut my DEMs, but I didn't know that an xarray tool allowed me to do that and also create circular ones! Thank you again!

[santisoler]

Hi @elisekazmierczak-geo! Sorry for my delayed reply. I cannot believe 3 weeks passed after my previous comment 😨.

Additionally, many points lie in valleys, and when I reduce the resolution of my DEM to discretize the topography into tesseroids, the interpolation places the surface above the actual terrain level. This resolution reduction is also dependent on the extent of my DEM, which brings me to my second question.

I'm afraid we are not doing anything special in those cases. I think what @craigmillernz suggested of using higher resolution grids around the stations is not only good to avoid such issue (observation point inside the source), but also to increase the accuracy of the terrain correction.

If I understood correctly, neither the prism nor tesseroid layer methods apply any special treatment to the area closest to the observation point and the prism-layer method doesn’t account for Earth curvature effects. I’ve tested several extension distances around the calculation point (and, consequently, different DEM extents), as the literature seems to vary quite a bit (e.g., 167 km, 22 km, etc.), and these choices can lead to differences on the order of several mGal.
Do you have a recommended extension distance specific to these methods?

You are right: the forward models for prisms and tesseroids just compute the forward model on the observation point, without any special treatment on areas closer to it. The prism forward is an analytic solution, so the result will be exact up to machine precision. On the other hand, there's no analytic solution for tesseroids, so we use a Gauss-Legendre Quadrature (GLQ) to approximate the integral with an adaptive discretization algorithm that splits the original tesseroids in areas closer to the observation point, so the quadrature is more accurate (you can read more about this algorithm in Uieda et al. (2016) and Soler et al. (2019)).

@craigmillernz technique for using lower resolution grids the farther away from the stations is great. I've wanting to implement an automatic version of that for a while now. I'll try to get back at it soon 🤞🏼

One of @leouieda's students will be working on terrain corrections. I think that will be also a nice opportunity to improve the tools we have in Harmonica for that.

[elisekazmierczak-geo]

Hi @santisoler , thank you for the time taken to answer me! 😊

Topographical correction is a tricky element because it is the only non-standardised correction and several methods exist. Even though Hammer zones are still widely used in literature, I was unsure about the relevance of their value in the Harmonica tools. Furthermore, this is the procedure that requires the most computation power and raises the question of which digital terrain model to choose. Our institute currently has 33,000 gravimetric measurement points spread over Belgium (30,500 square kilometres) which is a great opportunity, but requires a lot of computation power. That's why I first tried to make a prisms/tesseroids layer correction before testing local/global corrections. However, depending on the resolution and extension of the DTM, I have differences that can range up to 10 mGal per point in the more hilly areas of our country. I was hoping that there might be a more recent geostatistical analysis than Sprenke's (1989) for effective terrain correction, but I still feel that Hammer's values remain the benchmark for all types of terrain. I will therefore follow @craigmillernz 's advice 💪!

@santisoler , good luck implementing automatic topographic correction. It will certainly be a great step forward for your already very useful project!

Thank you both again for taking the time to respond, I hope to produce an as accurate as possible Complete Bouguer gravity anomaly map of Belgium! 🤓🇧🇪🍟