Significant differences in topography free gravity disturbance when extending the dem by 156km beyond the gravity data

[CampeolQ]

Description of the problem:

Hello, I'm currently using harmonica to process new gravimetric acquisitions in Belgium. I'm starting topographic correction using the prism method that you propose with harmonica. However, I'm wondering about the extent of the dem to be used in the prism_layer.gravity() method. In the bibliography, I read that it was not necessary to model prisms beyond 156 km around a station in this type of method, as their effect becomes negligible on the terrain correction. However, I found significant differences in the method's results when I extended my dem beyond 156km.

To illustrate my point, I've used your example of topographic corrections to gravity data from South Africa. I calculated the topography free gravity disturbance twice, integrating the 'padding' parameter. At this latitude, if we consider that 1° latitude = 100 km, I extended the prism grid by 156 and 175km around the data. I attached the notebook to this issue.

Normally, I'd expect the difference between the two results to be small, but strangely enough, it's sometimes of the order of twenty mgal. What's also strange from my point of view is the location of these differences, which are not on the edges. What do you think? Where could the error have come from?

Thank you in advance for reading and taking the time to reply.

Comparison_on_topographic_correction_with_different_buffer.zip

[craigmillernz]

Interesting. In some commercial terrain correction codes the code simply reflects the DEM edges out to the padding distance i.e. makes fake topography, so if your edge topography is extreme that might be continued to distance. Your dataset also extends offshore, so it might be an issue with bathymetry? (perhaps incorrect sign or something?). The highest change values also seem to be where the topo changes from "mountainous" to flatter?

Your two plots show the same area or am i mistaken?

Looking at your code it appears the padding comes from verde? (vd.pad_region), what DEM values are you applying to the padded coordinates?

[CampeolQ]

Thank you very much for your reply! I understand the idea of padding with fake topography, but in my case it's not this type of padding that's applied. As I mentioned above, I've mainly used the code in example Topographic Correction because in my own Belgian gravimetry project, I'm looking for the optimum extent of the DEM that I need to introduce for topographic corrections. I have therefore varied the extent of the DEM I introduce in the Harmonica example to see at what point the differences become negligible (in order to limit the DEM introduced and the calculation time).

I'm not sure this is a bathymetry problem. The example I've taken from Harmonica seems correct in this respect. In the code I've provided, I use Verde to create a zone of 1.5° or 1.75° extent, depending on the padding parameter, around the South Africa gravity data, and then cut out the chosen region from the South Africa DEM. If you compare the two plots, the altimeter scales are different and a peak at 3000m appears only on the right-hand plot.

Yes, indeed, the highest change values appears in the mountainous areas. Here is the complete DEM of South Africa, which I use to cut out the two areas of extent. The DEM extends 1.75833° east of the gravity data, so it's not a DEM limitation problem either.

[santisoler]

Hi @CampeolQ ! Sorry for the delayed reply. Thanks for opening this issue and for the detailed description of the problem you're facing. You also included a working example! Awesome!

After experimenting a bit with it, I think the source of the differences in the terrain correction is not coming from the forward modelling itself, but from the interpolation of the topography grid after its projection. Firstly, you define two different projections for each topography grid. These two projections are not the same because the lat_ts values are different (although they might be very similar). After that, each sliced topography is interpolated using the nearest neighbours method to get a regular grid of the topo heights in plain coordinates. The two grids end up having different values on points inside the survey area. These differences have an impact in the forward modelling of the prisms: since they are very close to the observation points, a difference in height can create a significant difference in the gravity field.

I modified your notebook to show you that the forward modelling is not the source of the issue. Instead of slicing and projecting the topography grid for each different padding, I projected the big topography grid (it might take a while to do it), and the sliced it for two different paddings. I ultimately computed the terrain effect of each sliced topography grid on the same set of observation points (which had been projected using the same projection used to project the topography).

If you take a look at the notebook I created, the differences between the two results are within the 0.04 mGal, which makes more sense than the differences you were seeing.

Comparison_on_topographic_correction_with_different_buffer_santisoler.ipynb.zip

I think this issue highlights the importance of correctly handling the projection and interpolation of the topography grid. One possible thing that I can think it could improve the results is to try the projection of the grid with another method, like the "linear". This is going to be much more expensive than the "nearest" (this is why we use the nearest in the example, we need it to run fast when we build the docs), but it might give more accurate results.

[CampeolQ]

Hello @santisoler,

Thank you very much for your reply and for the modified notebook you have attached. Indeed, I had not even considered the importance of correctly processing the projection before interpolating and padding the topographic grid. The scale's plot of the differences between the two results has become much more consistent. I'm also going to test interpolation using the "linear" method, thanks for the tips !

I think we can close this issue.