Pratt isostasy

[leouieda]

Description of the desired feature

PRs #17 and #36 added a function to calculate the Moho depth based on an Airy isostasy model. It would be great to have a Pratt isostasy calculation function as well. It would take the topography/bathymetry, respective densities, and a compensation depth and output the lateral density variation. The function should look something like this:

def isostasy_pratt(topography, density_crust=2800, density_water=1000,
                   compensation_depth=100e3):
    ...
    return density

The body will be very similar to isostasy_airy to allow numpy arrays or xarray.DataArray as input topographies. The docstring would benefit from having a figure to show what the terms in the equations are.

A good example case would be calculate the density variation along the mid-atlantic ridge.

[michiboo]

Hi Can I try tackle this issues?

[leouieda]

@michiboo of course, go right ahead 👍

[andieie]

Hi! I am attempting to solve this issue. To clarify, the output will be a grid of densities of ocean root and continent root? The input will be the assumed constant root thickness ? Does fatiando calculate an isostatic residual?
Sorry.. more questions than answers.

[santisoler]

Hi @andieie! I'm glad to have you onboard! I'll try to answer your questions.

To clarify, the output will be a grid of densities of ocean root and continent root?

Yes, but no necessarily a grid. The returned density array should have the same shape as the passed topography.
So, if topography is a 1d-array with scattered topographic heights, then density should also be 1d-array. If topography is a regular grid given as a 2d-array, then density should also be a 2d-array. But don't worry too much about the dimensions, Numpy will handle them automatically, so write it as topography is a 1d-array.

What @leouieda mentioned is that we should also enable the function to return an xarray.DataArray if topography is an xarray.DataArray, but don't be scared with this, you would basically need to copy some lines from harmonica.isostasy_airy().

The input will be the assumed constant root thickness ?

Yes, the function will assume a constant depth for the bottom boundary of the crust, as Pratt does. This depth should be passed through the compensation_depth argument.

Does fatiando calculate an isostatic residual?

Not specifically. I mean we don't plan to have a function that computes the isostatic residual anomaly entirely, but provide the tools to do so if anyone wants to. So, if I want to compute the gravitational field generated by the root of the crust (considering an Airy isostasy) or by the anomalous densities (Pratt isostasy), I would model it through prisms and then compute the gravitational field they generate (after #86 is merged, that will be very easy).

Sorry.. more questions than answers.

No problem! Feel free to ask anything! Hope my answers are helpful. If not, don't hesitate to write back.

You can follow Heiskanen & Moritz (1967) if you need some background for the math and/or a nice figure to understand what the terms in the equation are.