Building a real world model to calculate gravity

ASPECT
, ,
1

Hi all,

I would like to receive some input on how you would approach the problem I’m currently facing. I’ll start with a short description of what we are trying to do.

The Goal
The goal is to create a real world representation (3D) using S40RTS for the mantle and the winterc model for the top 410km. From this we want to use the gravity post processor and we want to do a Stokes solve.

Current Idea to Approach this Problem
Obviously, for the gravity calculation to be as good as possible we need the density to be in line with the density given in the winterc dataset. My best guess is to read in the density from winterc as ascii data in the Initial composition model. S40RTS can however, (as far as I know) only be read in using the Initial temperature model. Reading the winterc temperature dataset directly is also an option. However, the density obtained from this approach is not in line with the density given in the winterc dataset, so it will result in a different calculated gravity field.

Questions?

  1. How best to deal with both datasets to obtain realistic values for both temperature and density?

  2. Using the Initial composition model is only possible for the first 410km using the winterc ascii data, what to do with the rest of the mantle? Use the Initial temperature model?

  3. How to get from temperature perturbations (perturbed from PREM?) to absolute temperatures for the S40RTS data?

Hope my problem statement is clear enough, otherwise please ask for clarification. Thanks in advance!

Cheers,

Rens

2

Hi Rens,

My thought would be to use the ‘patch on S40RTS’ initial temperature condition. This allows you to use S40RTS but overwrite it down to a certain depth with an input ascii file. From your description, I think that’s what you want to do. The ascii read in has to be a shear wave perturbation so if you have a density perturbation field you can simply divide it by the “Vs to density scaling” that you set for your S40RTS model. That way, when you read it in and it multiplies it with the “Vs to density scaling” you get your density perturbation back out.
Your last question was about how you end up with absolute temperature. This depends on your material model so it’s a model choice you have to decide on. For the simple material models, the relationship is assumed to be linear so “ρ(T) = ρ_ref[1−α(T −T_ref)]” or “density_perturbation = (ρ(T) - ρ_ref)/ρ_ref= −α(T −T_ref)” where T_ref, rho_ref and alpha are parameters that you have to choose in your parameterfile (in the material model, see for example the simple material model).

I hope this is useful and helps you with your project!

Best,
Jacky

3

Hi,

Thanks for your reply! Do you know what kind of data the ‘Patch on S40RTS’ expects? Is it the Vs perturbation divided by the reference Vs values (d ln Vs), the absolute perturbation of Vs (d Vs) or the absolute Vs value? Furthermore, in case of one of the first two options: What is the reference Vs value? Is it PREM?

Thanks in advance!

Cheers,

Rens

4

Hi Rens,

I believe it’s the first, d ln Vs or (vs-vs_ref)/vs_ref. The reference vs depends on the tomography model, aspect doesn’t make explicit assumptions about that. It takes d ln vs and calculates d ln rho via a fixed or depth dependent scaling that you can prescribe in the input file. Absolute densities are then calculated from d ln rho using the reference density specified in the input file (material model).

Best,
Jacky

5

Hi Jacky,

Thanks again for your help!

I have another question about the way S40RTS is read in by ASPECT. Currently I’m looking at a python tool to read in spherical harmonics files and we are encountering some scaling differences in the order of a square root of two (plus some other small deviations). Are you a 100% sure the ASPECT implementation of S40RTS is correct? I saw in the /source/initial_temperature/S40RTS_perturbation.cc file that there are some notes about a square root of two scaling factor as well. We are not suggesting ASPECT is not working correctly. We just hope we can use the ASPECT implementation as a reference to benchmark our python code.

Cheers,

Rens

6

Hi Rens,

Apologies for this very late reply! S40RTS uses a different normalization than what (I believe) is normally used. Those are the comments that you will see in the source code and that might also be the reason why you get a different result when using an external package. As you point out, it’s related to a sqrt(2) prefactor. I originally implemented this incorrectly using a more standard normalization (sorry to be vague here, I don’t remember the details) and Shangxin Liu noticed shortly after that this doesn’t reproduce the perturbations shown in the paper and fixed the normalization. So I believe what is implemented is correct. A Sqrt(2) difference will be noticeable when you plot it so you can also compare your result and the aspect result to a figure in the original publication.

Sorry again that I’m replying so late and I hope that resolves your question but let me know if not! Also, any (potential) bug reports are always welcome :slight_smile:

Best,
Jacky