3D modeling and inversion of the density of a mountain

[shakasaki]

Description of the desired feature

Hello, I've just come across harmonica and I am still not very familiar with its use, so I wanted to ask whether what I am trying to do is feasible (or can be added as a feature that perhaps I can also support).

Basically, I have a 3D DEM model of a mountain. There is also a tunnel going through this mountain, and I will take gravity measurements along this tunnel. I want to use these for an initial inversion, which I will then use to design a more elaborate survey of taking measurements on the mountain itself.

Ideally, I would want to create a 3D grid that uses the DEM of the mountain and the tunnel geometry (see figure attached, which is obtained using pyvista, of the geometry I am talking about).

If I understand correctly, the way to do this currently is to use the harmonica.prism_gravity function.
I saw thay pyGimli is supported. IS it possible to create the mesh using pyGimli?
Or is there another meshing tool supported in harmonica?

And how do I then run the inversion? Apologies for the many questions but I could not find this info in the documentation!

Thanks for the help, and if there is some functionality that I can help with, I would be happy to do so!

Extra info: Just to get an idea of the scale, the DEM model spans roughly 6x6km of surface, the highest peak of the mountain is 3.2km and the tunnel is 5km long

Are you willing to help implement and maintain this feature? Yes

[welcome]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible.

You might also want to take a look at our Contributing Guide and Code of Conduct.

[santisoler]

Hi @shakasaki! Glad to see that Harmonica is getting more and more attention.

On its current development stage, Harmonica doesn't have any inversion algorithm, although we do have the plan for future implementations of the types of inversions that the old fatiando package has. Nevertheless, we are more interested in non-conventional inversion algorithms. There are some very well designed packages for 3D meshes inversions, like SimPEG and pyGIMLi, so we won't intend to overlap any Harmonica feature with any of those. If you want to read more about the design decisions that Fatiando took a couple of years ago, please read this blog post.

On the other hand, we do have functions for computing the gravity effect of prisms and we are actually working on a function to create a layer of prisms (see #186), which top and bottom boundaries can be modified by the user to (for example) drape the topography of a DEM or the geometry of the Moho. In conclusion, we are working on simple ways of computing the gravitational effect of a given terrain by approximating it through prisms, but not for carrying out an inversion.

Just for curiosity, your intention is to define a mesh of prisms and invert their density?

[shakasaki]

Hi @santisoler and thank you for the quick and detailed response!

Great that you pointed out the existing projects. I was aware of pyGIMLi but not SimPEG. And while the former only has 2D gravity modeling, the latter seems to do it also in 3D. I'll definitely check it out!

My intention is to have a forward solver that can handle both topography and also data within the mesh.
That is, I am collecting gravity measurements on a mountain-top but also within a tunnel that crosses the mountain. So, discretizing densely close to the tunnel, and on the measurement points on the surface, would be ideal.

I suppose compared to the SimPEG forward solver, one advantage of Harmonica is the use of prisms instead of the regular grid type (even though the SimPEG grid seems to be refinable). Then I could also use the exact tunnel surface (it's a bit rugged) that I have from a laser point cloud.

Setting up the inverse problem to get the underlying density is the final aim here. But even if that is not currently present in Harmonica, it would be a part that I can contribute to.

[santisoler]

Closing this issue. Feel free to comment if you have any update.