3D gravity inversion of prism layers

[mdtanker]

Description of the desired feature:

I'm working on a 3D gravity inversion for my Ph.D. to model the bathymetry beneath a float ice shelf. I'm using a vertical-prisms approach, and I'm mostly interested in a geometry inversion, where each prism's density remains unchanged, but the prisms' tops or bottoms are updated to minimize the observed gravity - forward gravity misfit. While this inversion is in the early stages, I'm already using several harmonica tools so I think it would be great to include a fully supported inversion module in harmonica.

Below is a description of my specific use case, a list of what I already have coded, and ideas for additional features to add.

Current state of inversion

Model setup

• input an arbitrary number of layer topographies (i.e. ice surface, ice base, bathymetry, basement, Moho)

  • xarray.load_dataset() to load netcdf's as grids
  • pygmt.grdsample() to make sure all grids' extents match

• create a layer of vertical prisms between each set of grids

  • harmonica.prism_layer()
    • example: for water layer: surface = ice base, reference = bathymetry
  • assign a density to each prism (currently using constant density within each layer)
    • prism_layer(... properties={'density':water_density)
  • allow different cell sizes between layers
    • if lower grid cell size doesn't match upper grid's, use pygmt.grdtrack() to sample the lower grids values at the grid's prism locations

• calculate forward gravity of each layer

  • harmonica.prism_layer.gravity() for each prism layer

• calculate gravity misfit

  • add forward gravities for all layers, and subtract from observed gravity

run geometry inversion

• designate an active_layer ex: active_layer='bathymetry_prisms'
• each iteration of the inversion will yield a surface_correction
  • set a max surface change per iteration
  • applied to the active_layer tops, as well as the bottoms of the layer above active_layer
    • harmonica.DatasetAccessorPrismLayer.update_top_bottom()
  • recalculate the forward gravity of these 2 updated layers.
    • recalculte the misfit with these updated forward gravities before starting the 2nd iteration
• currently this inversion works by finding the least-squares solution to the matrix equation Ax=b, where A is the Jacobian matrix of the vertical derivative of gravity acceleration, and b is the initial misfit between observed and forward gravity.
  • the least squares solution is found with scipy.sparse.linalg.lsqr()
  • the vertical derivative is approximated with the Hammer approximation of prisms, using annuluses as opposed to prisms.

This is just the method I already had implemented, suggest from my advisor, but I'm happy to alter this and would appreciate suggestions.

Features to add:

• enforce constraint point (points of known bathymetry) by multiplying surface_correction by a constraints grid * 0's at constraints, linearly increasing to 1's at a specified distance from nearest constraints * scipy.spatial.cKDTree.query() to get grid of minimum distances to constraints * gmt grdmath LE and GE to set constraints to 0 and other points to 1 * gmt grdblend to merge the clipped grids * and rioxarray.interpolate_na() to linearly interpolate between 0 and 1

• allow irregular grid spacing within each layer (i.e. grid spacing increase for a buffer zone, outside the main area of interest)

  • could use discretize.TreeMesh(). from SimPEG.

• allow depth-dependent density within each prism. Only for the forward modelling

  • would be useful for firn/snow compaction into ice, and sedimentary basins.
  • similar to https://github.com/fatiando/harmonica/pull/269

• density inversion

  • similar to the above geometry inversion, but updated each prism's density to minimize the misfit.
  • this is useful for the initial setup of the geometry inversion, allow the regional field to be accounted for in a layers density, leaving only the residual field left to invert on.

Discussion

• Do we want to refer to this type of inversion as a geometry inversion, structural inversion, boundary inversion or something else? There doesn't seem to be much of a consensus on what you call this.
• I'm pretty new to inversion theory, so if anyone has recommendations on where to learn this material please let me know!
• I think some of the packages I have linked above could be replaced with Fatianado packages to keep dependencies minimal.
• All the inputs in my inversion are stored in a nested dictionary layers, with a dictionary for each layer, which includes a .nc filename, a resolution, constant density values, and a dataframe and dataset for each layer. This makes it easy to apply functions to an arbitrary number of input layers (example below). If there's a better method for this let me know.

layers_list =[ 'ice', 'water', 'bathymetry', 'basement',]
spacing_list = [10e3, 10e3, 10e3, 10e3,]
rho_list = [ 920, 1030, 2600, 2800,]
fname_list=[
            '../inversion_layers/BedMachine_surface_filled_inv.nc',
            '../inversion_layers/BedMachine_icebase_filled_inv.nc',
            ../inversion_layers/BedMachine_bed_inv.nc',
            '../inversion_layers/ROSETTA_basement_BedMachine_bed_inv.nc',
            ]

# create nested dictionary of layer properties 
layers = {j:{'spacing':spacing_list[i], 
            'fname':fname_list[i], 
            'rho':rho_list[i]} for i, j in enumerate(layers_list)}

for i, j in enumerate(layers_list):
    # load each .nc as a xarray dataset
    layers[j]['grid']=xr.load_dataset(layers[j]['fname'])
    # create prism layer between each set of input grids
    layers[j]['prisms']=hm.prism_layer(
        coordinates=(list(layers[j]['grid'].x), list(layers[j]['grid'].y)),   
        surface=layers[j]['grid'].z, 
        reference=layers[reversed_layers_list[i-1]]['grid'].z,
        properties={'density':layers[j]['grid'].density})

# calculate forward gravity for each prism layer
for k, v in layers.items():
    df_grav[f'{k}_forward_grav'] = v['prisms'].prism_layer.gravity(
        coordinates = (df_grav.x, df_grav.y, df_grav.z),
        field = 'g_z', progressbar=True)

• to remove edge effects, I'm using a buffer, currently 200km in each direction, as shown in the below figures by the black square. You can see the edge effects in the second to last panel, which is the lowest layer of prisms. I am doing all calculations on the full region, but only showing results within the buffer zone. Does this seem like a reasonable approach? It adds a lot of prisms, but this could be mitigated by using discretize.TreeMesh(). for the buffer region. I'm unsure if I should be calculating misfit for the whole region or just the inside. Also, I'm not sure if I should invert the prisms just within the inside, or for the whole region.
  

Are you willing to help implement and maintain this feature?

Yes I'm happy and excited to implement this, but I'm in the last year of my Ph.D. so the outcome of this inversion is more pressing than the implementation of the inversion into Harmonica. Hopefully they can happen simultaneously. @santisoler has suggested looking at both the old Fatiando module fatiando.inversion and some code related to Uieda and Barbosa 2016 for inspiration on how to implement this, which I will do.

[LL-Geo]

Just find this paper, that might give an idea of how to include different density contrast (maybe different geology) into the inversion. Haas et al 2020. https://academic.oup.com/gji/article/221/3/1896/5808816?login=true

[craigmillernz]

This would be very useful feature. There is some overlap with SimPEG here. Perhaps SimPEG inversion engine could be used and linked to Harmonica for forward model with prisms. I think @santisoler is looking at this in SimPEG at the moment to rework SimPEG gravity forward modeling. It would save reinventing a lot of inverse code that is already well used and supported.

VPMG commercial code that does this style calls this a geometry inversion, VPMG also has the option to alternate 1 geometry iteration and one property (i.e density) iteration then back to geometry iteration etc to solve for both geometry and density.

[santisoler]

Hi @craigmillernz. Thanks for joining the conversation. I agree that we should make efforts to avoid overlaps with SimPEG. And you're right, I'm working on building some common ground for prism forward modelling between SimPEG and Harmonica. I'm looking forward to getting more awesome results.

Nevertheless, I like the way that a geometry inversion could be built on top of Harmonica. SimPEG is mainly developed for inverting the physical properties of an already discretized space. I'm guessing that trying to implement something like the geometry inversion carried out by Uieda et al (2017) on top of the SimPEG framework could be challenging. But again, I don't mean to write a complete separate inversion framework in Harmonica just for that. Maybe we could start with specific inversion functions, like invert_moho for inverting the Moho or invert_basin for a sedimentary basin, more in sync with a use-case oriented development. What do you think?

[mdtanker]

Thanks for pointing that out @LL-Geo, I'll look through it.

[mdtanker]

Santi I like your idea of having a few specific inversions for common tasks, but it might end up with be a lot of duplicated code. What do you think are the benefits over having a generic invert_prism_geometry function, with specific gallery examples for applying it to common use cases?

[craigmillernz]

I tend to agree with Matt that a more generic function with use cases might be better. I was really just wondering how easy it would be to couple the prisms discretization onto the simpeg inversion framework, You'll know best about that! If they aren't really compatible, then a separate inversion within Harmonica would be appropriate.

[craigmillernz]

Matt you might find this useful for general inversion theory - from UBC. This used to be a standalone website, but it looks like you now have to download. Follow the instructions on the page. https://gif.eos.ubc.ca/IAG

[mdtanker]

Hi everyone, sorry for the lack of updates here! I've been struggling to implement regularization for a while and was wondering if anyone had ideas to help. I've added some smoothness regularization within the least-squares solver with a damping term (with verde.base.least_squares), which helps stabilize the inversion. But, I'd like to add some hard constraints into the regularization as well, which seems to be referred to as "smallness" regularization (as discussed here).

There are a series of seismic constraints on the bathymetry depths which I would like to be unchanged (within error margins) after the inversion. Ideally, at each iteration, the correction value at these points would be 0, and the smoothness term will help ensure these constraint points don't turn into pedestals, as the unconstrained bathymetry surrounding them changes.

Does anyone have ideas for how to implement this? Would it be within the least squares solver?

I've tried to follow along with similar implementations in PyGIMLi, SimPEG, and the legacy Fatiando, but I've either only seen smoothness regularization, or I'm unable to understand how it's coded.

Any help or suggestions would be appreciated!

[santisoler]

Hi @mdtanker. I'm glad you are getting into the rabbit hole of inversion, it's a fun one haha!

The good news here is that you are already using some simple kind of smallness regularization. The verde.base.least_squares function relies on the sklearn.linear_model.Ridge class that minimizes the objective function:

$$ \phi(\mathbf{m}) = \lVert \mathbf{d}_\text{obs} - \mathbf{G} \mathbf{m} \rVert^2 + \alpha \lVert \mathbf{m} \rVert^2 $$

The second term that involves the alpha (damping parameter) and the l2 norm of the parameters is a Tikhonov zeroth order regularization. SimPEG implements the smallness by defining this regularization term as the l2 norm of the difference between the model parameters and the reference model. In this reference model you can include you constrains:

$$ \phi(\mathbf{m}) = \lVert \mathbf{d}\text{obs} - \mathbf{G} \mathbf{m} \rVert^2 + \alpha \lVert \mathbf{m} - \mathbf{m}\text{ref} \rVert^2 $$

Here's a nice chapter to get into these topics: https://pubs.geoscienceworld.org/books/book/2111/chapter-abstract/114879875/Inversion-for-Applied-Geophysics-A-Tutorial?redirectedFrom=fulltext

The problem is that you won't be able to use the Ridge class for defining a reference model and you'll probably have to write you own implementation of the regularization.

There are also some lecture notes from @leouieda and @birocoles. They have really nice stuff there, particularly if you need to implement this stuff: https://www.semanticscholar.org/paper/T%C3%B3picos-de-invers%C3%A3o-em-geof%C3%ADsica-OliveiraVanderlei-Leonardo/fd651157443f37ad603057a10d683322ff123263

[WYJLCYWHZ]

Hi @mdtanker ,I came across your project on GitHub and I'm truly impressed by the remarkable feature you've been working on. Your research on 3D gravity inversion for modeling the bathymetry beneath a float ice shelf sounds fascinating. I'm particularly intrigued by your implementation of the vertical-prisms approach and the focus on geometry inversion to minimize the observed gravity misfit.

I'm genuinely interested in learning from your implementation and gaining insights into the techniques you've employed. It would greatly benefit my own learning and allow me to explore potential areas for improvement.I would be grateful if you could kindly share the code you've already developed. Thank you (^▽^)

[mdtanker]

Hi @WYJLCYWHZ. apologies for the slow reply. Thank you for your kind words on my work. The inversion code is included in my repository RIS_gravity_inversion, but it is written specifically for my use-case and is not very well documented or tested yet. In the coming months, I will work on making this inversion code more accessible and understandable for other users. I'll try and post updates in this thread or in a Fatiando discussion. Let me know if you have specific questions.

[WYJLCYWHZ]

Hi @WYJLCYWHZ. apologies for the slow reply. Thank you for your kind words on my work. The inversion code is included in my repository RIS_gravity_inversion, but it is written specifically for my use-case and is not very well documented or tested yet. In the coming months, I will work on making this inversion code more accessible and understandable for other users. I'll try and post updates in this thread or in a Fatiando discussion. Let me know if you have specific questions.

Thank you very much (^▽^)