unrealistic temperature field at slab tip during subduction

[ythaha]

Hi there, I'm a beginner in using underworld and i found it nice to use. But i encounter a problem. I'm reproducing Holt 2015 , a simple subduction model : using weak crust ( thickness of 15 Km) to decouple slab and overriding plate, all boundary are insulated and free slip. Initially slab and both (subducting and overriding) plates have uniform temperature 0 (non-dimensional), other region (mantle) has uniform temperature 1.

The initial condition seems ok, both temperature and viscosity field.

although the slab tip is somehow wrong —— viscosity is larger at the tip

My problem is, during subduction, temperature at slab tip get unrealistic —— with higher (1.8) and lower(-1.8) temperature at the slab tip. This caused much larger viscosity and finally it crashed, maybe the viscosity is too much large. ( results are shown at comment below)

I think the reason is there is some wrong setting caused by misunderstanding the usage of underworld. Below is my script.

Thank you for any instruction!

https://gist.github.com/ythaha/116d785afe97055f97549ec87951e1a2

[ythaha]

sorry i forget the temperature field. here it is, field at 150 step:

and viscosity

[julesghub]

Hi @ythaha, Welcome to the Underworld! Great work setting up this model, I'm unfamiliar with it but it appears it should be reproducible.

I reckon the "oscillations" you see at step 150 are due to advecting the gradient (0, 1) at the tip of the slab. Underworld's current AdvectionDiffusion numerical method (SUPG) can produce these oscillation when large field gradient's exist in an advection dominate model.

A workaround is to allow the temperature field to diffuse somewhat before allowing the model (i.e. materials+temperature) to advect. The aim would be to smooth the temperature gradients over an element or two as opposed to a sharp change in value between single elements.

To do this you can: 1) zero the velocity field: velocityField.data[...] = 0.0 2) Pick a dt to diffuse the temperature gradients and run only the advDiff system. ie

dt_diffuse = uw.scaling.non_dimensionalise((xxx * u.year)
advDiff.integrate(dt_diffuse)

with the zeroed velocity field. 3) output the temperature field and inspect the suitability of this diffused starting conditions.

Give that a try and report back.

[ythaha]

Thank you @julesghub for your guide! And I did some tests following your instruction.

1. first, i tried to diffuse the temperature field for 1Ma/10Ma before advect material and temperature, it does look better. But still had some oscillations at later timesteps.
2. then, i tapper the weak crust, limiting it to the top 150 Km (as many papers did) . The oscillations disappear, but a cold region formed at the tip of slab, which i think is unrealistic

So i guess the oscillation might be caused be the low viscosity layer on top of slab, although i don't know why.

3. finally, as you mentioned SUPG, i test using SLCN for AdvectionDiffusion solver, the result seems better, and the timestep between each step is larger, which means the modeling is faster.

Indeed i'm afraid that the timestep maybe too large. Using SLCN, the timestep is about 2 times larger than SUPG. (SLCN timestep is about 45K year on average)

So, can i use SLCN for future modeling? and do i need some particular settings for this method to get right results? (such as timestep or so.)

It would be a great help if anyone can provide some information about this method. Thanks!

[ythaha]

do i need compare max timestep for swarm advection and temperature? like this in UWGeodynamics

self._dt = 2.0 * rcParams["CFL"] * self.swarm_advector.get_max_dt()

if self.temperature:
    # Only get a condition if using SUPG
    if rcParams["advection.diffusion.method"] == "SUPG":
        supg_dt = self._advdiffSystem.get_max_dt()
        supg_dt *= 2.0 * rcParams["CFL"]
        self._dt = min(self._dt, supg_dt)

https://github.com/underworldcode/UWGeodynamics/blob/e99e225d4b149497827e09d14b3f4df5822b380a/UWGeodynamics/_model.py#L1603

[ythaha]

do i need compare max timestep for swarm advection and temperature? like this in UWGeodynamics

self._dt = 2.0 * rcParams["CFL"] * self.swarm_advector.get_max_dt()

if self.temperature:
    # Only get a condition if using SUPG
    if rcParams["advection.diffusion.method"] == "SUPG":
        supg_dt = self._advdiffSystem.get_max_dt()
        supg_dt *= 2.0 * rcParams["CFL"]
        self._dt = min(self._dt, supg_dt)

https://github.com/underworldcode/UWGeodynamics/blob/e99e225d4b149497827e09d14b3f4df5822b380a/UWGeodynamics/_model.py#L1603

it seems I should do this to get right max dt

[ythaha]

2. then, i tapper the weak crust, limiting it to the top 150 Km (as many papers did) . The oscillations disappear

with SLCN, even without limiting weak crust depth, the temperature still had no oscillations and looks well.

As I know, citcom/citcomS also use SUPG for advection, but I didn’t encountered something like this using them, or maybe I had some misunderstandings about that?

[jmansour]

Hi @ythaha

With regards to the timestep size, note that the SUPG timestep is a necessary maximum for stability. Indeed, you should find that if you decrease the timestep size, the oscillations should be reduced. I can't comment on citcom/citcomS, but possibly it defaults to a more conservative timestep size, although there are other choices for the SUPG implementation that can have stability implications.

The SLCN method is however unconditionally stable, and you are free to choose whichever timestep size you wish, but note that something around the CFL determined size is usually best for accuracy. This is the same for the swarm advection, which is unconditionally stable, but again for accuracy should be around the CFL size.

[ythaha]

Hi @jmansour , thank you for your information about CFL. But here i have another problem:

I tried using SLCN for modeling. At first, it looks well. But After adding phase change function, the temperature field show some "peaks".

This is temperature with phase change function

temperature along a horizontal line at 1000 km depth. scaling is 1400K.

From this pattern, i suspect the peak corresponds to local boundary. Since the i use 28 cores, the mesh maybe split to 14x2 local meshes. ↓

    Global element size: 1024x256
    Local offset of rank 0: 0x0
    Local range of rank 0: 74x128

This is without phase change, which show no peaks.

---

Here is my code for phase change. (The full code is here : https://gist.github.com/ythaha/51008a41295e252d34946552b8211f98 )

gamma_410 = 3.0e6 * 1400 / 3300 / 10 / 2e6
temp_410 = 1.07462
delta_rho_410 = scaling.non_dimensionalise(273. * u.kilogram / u.meter ** 3)
width_410 = scaling.non_dimensionalise(13 * u.kilometer)
Phase_change_410_func = -1. * fn.input()[1] - 410/2e3 - gamma_410 * (tempField - temp_410)
Phase_change_410_density = delta_rho_410 * 0.5 * ( 1. + fn.math.tanh(Phase_change_410_func/width_410))

# %%
force_fn = (
    0.0,
    (uw.scaling.non_dimensionalise(refDensity)
    * (1.0 - uw.scaling.non_dimensionalise(3e-5 / u.degK) * tempField)
    + Phase_change_410_density 
    )
    * uw.scaling.non_dimensionalise(-10.0 * u.meter / u.second ** 2),
)

The formulation is from this paper

---

Then i test SUPG, using 0.6 * get_max_dt

    dt1 = advDiff.get_max_dt()
    dt2 = advector.get_max_dt()
    dt = 0.6 * np.min([dt1,dt2])

The temperature field seems normal.

(although upper mantle show too vigorous convection, which i will check later)

---

Have you ever experienced something like this unusual "peak" temperature? I'm happy with SLCN before i encountered this problem, so i'd like to learn more about this. Thanks for any help you can provide.

[julesghub]

@ythaha, I agree with you, the "peaks" do look like parallel issues. Is the mesh deformed? This may explains the results as the SLCN algorithm isn't implemented for arbitrary deformations to the cartesian mesh. For arbitrarily deformed meshes we recommend SUPG.

[lmoresi]

I would be a little bit careful here as the SLCN approach “reaches back” through the temperature field to the launch point of the point that is being computed. I can see that introducing some complexity into whether the phase change is considered if the path crosses the phase change

I would also think that it would be necessary to make sure the shadow zones reflect the correct temperature. That would introduce some parallel issues that could be different between the two methods.

L

Prof Louis Moresi @.*** www.moresi.infohttp://www.moresi.info/ www.geo-down-under.org.auhttps://www.geo-down-under.org.au/ @LouisMoresihttps://twitter.com/LouisMoresi

On Wed, Jul 7, 2021 at 10:24 AM, Julian Giordani @.***> wrote:

@ythahahttps://github.com/ythaha, I agree with you, the "peaks" do look like parallel issues. Is the mesh deformed? This may explains the results as the SLCN algorithm isn't implemented for arbitrary deformations to the cartesian mesh. For arbitrarily deformed meshes we recommend SUPG.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/underworldcode/underworld2/issues/544#issuecomment-875171444, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADABPI3TARI4TAHDMMX5DX3TWONDPANCNFSM45VQA5JA.

[jmansour]

@ythaha can you share your model so that we can reproduce the issue? Attach it here, or email, or whichever mechanism you're comfortable with.

[ythaha]

Yes, here is my script that produced those peak. test_5.py.txt ( before running, please create the corresponding output directory)

Recently i found that those peak well be suppressed a lot by adding adiabatic heating in model, which will also suppress small convections at upper mantle. Those convections at upper mantle maybe caused be density instability by having mantle temperature following adiabat ( temperature increase with depth ).

Here is my code with adiabatic heating. i'm not confident with my implement, it will be a great help for me if you can check it. Thanks for all the help from your team! test2_1.py.txt

[ythaha]

@ythaha, I agree with you, the "peaks" do look like parallel issues. Is the mesh deformed? This may explains the results as the SLCN algorithm isn't implemented for arbitrary deformations to the cartesian mesh. For arbitrarily deformed meshes we recommend SUPG.

Yes, the mesh is deformed to increase resolution around subduction zone and upper mantle.

model_length = 1e4 * u.kilometer
model_top = 0.0 * u.kilometer
model_bottom = -2e3 * u.kilometer

mesh = uw.mesh.FeMesh_Cartesian(
    elementRes=(1024, 256),
    minCoord=(-0.5 * model_length, model_bottom),
    maxCoord=(0.5 * model_length, model_top),
)

with mesh.deform_mesh():
    mesh.data[:, 0] = mesh.data[:, 0] * np.exp(
        (np.abs(mesh.data[:, 0]) - 0.5 * model_length) / 3.0
    )
    mesh.data[:, 1] = mesh.data[:, 1] * np.exp(
        (np.abs(mesh.data[:, 1]) - np.abs(model_bottom)) / 1.5
    )

[julesghub]

https://github.com/underworld-community/ythaha_oscillations

@ythaha and @jmansour here is a repository for working on the model.