Suggestion: seismic velocity to density conversion

[cmeessen]

Dear all,

for my research I am interested in how we can constraining mantle densities from seismic tomography. Therefore, I implemented a python version of the conversion approach by Goes et al. (2000) who proposed a way to convert s/p-wave velocities to density and temperature. The authors use Newton iteration to match synthetic and observed velocities which makes a direct implementation in Python very slow. Since I was applying the conversion on 3D regular meshes, I was able to accelerate the conversion by using lookup-tables for the individual depth values of the mesh. With about 10 000 velocity values and 8 depth slices it runs in less than 10s. You can have a look at the code at this repository, I hope the comments make the code understandable ;).

Since fatiando is using PrismMeshes, I thought this conversion might be a nice addition to the package. @leouieda already suggested to create a new petrophysics module or similar.

So far, this would be a rough roadmap of what needs to be done and what could be improved:

• update syntax from Python2 to Python3
• implement function to read properties from PrismMeshes
• implement function to define mantle rock mineralogy (at the moment uses a csv file), maybe with predefined mineralogies
• implement function to write densities and temperatures back to the PrismMesh
• clear out unused functions (or keep a command line interface?)
• documentation
• tests
• (improvement) alternative pressure calculation with PREM or other reference models
• (improvement) parallelise for large data sets
• (improvement) allow for compositional heterogeneity by storing the mineralogy in the PrismMesh. This would be kind of a very distant stretch goal and should be restricted to a limited number of different mineralogies (e.g. Harzburgite and Basalt for mantle and subducting slab), otherwise there would be the need to calculate too many lookup tables

Except the parallelisation and compositional heteorgeneity these tasks are not too difficult. So what do you think? I'd be happy to answer any questions.

Looking forward to your comments!

Cheers,
Christian

[leouieda]

@cmeessen thanks for sharing your code!

I think a new fatiando.petrophysics module would be a good place to house physical property conversions. I'm open to suggestions for alternatives.

If we're going with a new package, it would be better to create the empty package in a separate PR first to make things easier. It's not a lot of work but the package needs to be included in a few places in the docs. I can do this quickly.

When porting your code, it's best to do it in small portions so that we can quickly test, document, and merge. I have a tendency to start big PRs and that causes a lot of trouble when the work ends up taking a long time. I always have to split it up later. So don't do what I do :)

Regarding the command-line, I think that is best kept out of Fatiando. It could be a separate program that uses Fatiando to implement an application. We can host it on the organization account if you want and have it linked in the site. But that will give you more flexibility with the releases and updates.

cc @fatiando/contributors

[victortxa]

I haven't looked the code of @cmeessen, but something I was wondering sometime is if Bruges could be used in Fatiando (I don't know how Bruges is going, but I know there are Petrophysics and Rock Physics modules there). I'm mentioning this because of the possibility of starting this Petrophysical module in Fatiando (should we use what is already done? How?).
What do you think?

[leouieda]

@victortxa that's a good point. I don't know how active bruges is or their status. I haven't seen updates for a while. Cc @kwinkunks @ben-bougher @EvanBianco

That's the problem with trying to have monolithic packages like Fatiando. I get their initial appeal but it gets frustrating when you want to use one part of the library but not all the rest that comes attached to it. Might be time we consider splitting Fatiando into several packages. They could all live under the Fatiando umbrella still.

[EvanBianco]

@victortxa, Bruges hasn't had much added to it's compilation for a while, but I think it might be the right place for Petrophysical equations and calculatory stuff. It's heavier in rock physics and seismic equations so far, but that hasn't been deliberate. We've just started with the things we know best and are trying to fill it out from there. Happy to discuss / map out a Petrophysical sub-module within Bruges and we can start filling in the gap we have there.

[victortxa]

@leouieda any thoughts about this (integrate the Petrophysical package from Bruges)? Maybe a general guide for this kind of thing can help in anothers future modules.