Effective viscosity in the Viscoelastic formulation

ASPECT
1

Hello everyone,

I have been working on the effect of elasticity in RT type drips. To do this I have run some simple tests checking the results of the models through the topography and the maximum value of the second invariant of the deviatoric stress tensor \sigma'^{2nd}. However, the results are not in agreement with published results.

I tried calculating \sigma'^{2nd} using the full deviatoric stress tensor (2\eta\dot{\epsilon_{ij}}+\sigma_{ij_{elastic}}) using its components to calculate the invariant, and I also calculated \sigma'^{2nd} from the principal stresses using the recent PR for their visualisation, but both results are off. Also, the topography which comes directly from ASPECT’s solution tends to produce more surface deflections than the expected.

I am wondering if this difference is produced by the formulation of the effective viscosity: \eta_{eff}=\eta\frac{\Delta t^{e}}{\Delta t^{e}+\frac{\eta}{\mu}} , which is constant throughout the entire model while for timescales comparable or grater than the relaxation time \frac{\eta}{\mu}, the effective viscosity should approach the “true” viscosity with a formulation such as: \eta_{eff}=\eta-(\eta-\eta_{eff_{0}})e^{-\frac{\mu}{\eta}t}.

I would like to know if this is a simple thing to test maybe by modifying the equation for \eta_{eff} in a material model plugin or if there are more difficulties.

Thanks a lot!
Daivid

2

Hi David,

Thank you for posting to the forum! @bobmyhill and I discussed your question earlier today, but could not think of a reason why the existing viscoelastic implementation would not be giving the correct stress output.

As such, we think the difference between the ASPECT results and the ones you are comparing with are likely due to the different effective viscosity formulation.

It would not be hard to compute a different effective viscosity based on the formula you posted. Shall we have an offline discussion about how to add that in? We can also discuss that the next ASPECT user meeting in ~ 2 weeks.

For reference, can you add a link to the specific study you are comparing results with? I recall it was work by Boris Kaus, but wanted to double check. Also, just to be clear, is it correct that the results you are comparing with are not from an analytical solution?

Thanks!
John

3

Hi John,

Thanks a lot for your reply, it would be great to discuss how to implement this, please let me know when could we talk about it! I will also attend to the next ASPECT user meeting.

The link to the work by Boris Kaus can be found in this link. The results I am comparing with are his numerical modelling results (See figures 9, 10 and 11).

Regards,
David

4

Hi David,

Thanks for the link. Yes, I think the difference in the effective viscosity formulation is likely the key difference.

Let’s discuss how to implement this at the next ASPECT user meeting, as others in the community may be interested as well.

Cheers,
John

5

Hi David,

Not much to add to John’s reply, but here are a couple of notes related to the last paragraph of your reply which may help stimulate further discussion:

  • The “effective viscosity” provided by the function calculate_viscoelastic_viscosity in elasticity.cc is not \eta in the usual viscous expression \sigma_{ij}=2 \eta \dot{\varepsilon}^T_{ij}. Indeed, there can be no such expression for viscoelastic materials. A good illustration of why this must be true: viscoelastic relaxation of stresses can take place when the total strain rate is exactly zero (in which case this constitutive form would not be useful). Even if strain rate is non-zero, \eta will generally be of fourth rank: \sigma_{ij} \propto \eta_{ijkl} \dot{\varepsilon}^T_{kl}, and of second rank even for isotropic viscoelastic materials.
  • The calculate_viscoelastic_viscosity function actually returns the effective viscosity for the very special case that the material was unstressed at the beginning of the timestep. This is not generally the correct viscosity, so some corrective terms are required.
  • As I understand it (not having worked with the elastic rheology yet), those corrections are supplied by the elastic and reaction outputs in elasticity.cc. These are used by other parts of the ASPECT code to allow the preexisting elastic forces to be added to the governing equations.
  • The ASPECT implementation follows the Moresi et al (2003) paper cited in the code.

I hope this helps a little bit.
Bob

6

Hello everyone,

I have done a few modifications to the code to include the desired time dependence in the effective elastic viscosity. When I run the code it produces the desired change and it runs fine for the first few time steps. However, at some point the time step starts decreasing finally approaching 0. Hence the model stalls (It does not crash, instead it reaches a stagnation point where it tries to progress but the time step is very close to 0). I have been trying to understand the cause of this but I have not been able to find the problem.

I could submit a PR with the modifications that I have done so that you can check it, but since it doesn’t work yet and I am not sure if this change in the formulation of the effective viscosity will work, I don’t have all the tests/benchmarks and documentation that you ask for in the “PR checklist”.

Do you know what could be causing the problem with the time step? Should I submit the PR or show you the modifications to the code in another way?

Thanks a lot!
David

7

Hi David,

Thanks for posting!

My first guess would be that a very large decrease in viscosity (or large increase in elastic forces) is leading to very small time steps. Can you post a link to the following?:

  1. Your GitHub branch with the changes?
  2. A folder (google drive, dropbox, etc) with the model results? Alternatively, feel free to post some images here.

Thanks!

Cheers,
John

8

Hello David and John,
I don’t know the exact problem that you are working on, but I’ve been building a visco-elastic subduction model and choosing an appropriate fixed elastic time step for your particular problem is a little bit tricky. A smaller value leads to a lower effective viscosity (there’s a nice note in the code about this: eta-elastic = t_e*mu, the effective viscosity is given by 1/eta_eff = 1/eta-elastic + 1/eta-visc. In setting up the subduction model, I first chose a t_e of 1e5 (100 kyr) with rigidity of 75 GPa, this reduce my maximum viscosity of 1e24 Pa-s to something like 5e23 Pas and given the driving forces from the slab, this meant I had plate velocities of up to 1-5 meters per year (yikes!), which requires a tiny time-step. Setting a fixed numerical time-step helped this somewhat, but not enough. Increasing t_e to 1e6 yrs reduced to 30 cm/yr and now I’m also adjusting a few other things. (this time-scale is appropriate as I’m interested in the deformation on the time-scale of material moving through a deforming region in the slab of about 100-200 km lengthscale). I’m also using a fixed numerical time-step that is much smaller that t_e (only up to 2000 years, for testing, but I also had problems letting the CFL condition set the numerical time-step. The other thing I found is that inconsistencies in the model setup lead to elastic stress anomalies (view the stress_xx, stress_yy, stress_xy fields to see these). For example, at one point I had these along a crust/lithosphere compositional/viscosity boundary that needed to be better resolved. Maybe this helps. - Magali

9

Hi John,

Thanks a lot for your reply! In my case I am producing a fast increase in the viscosity, could this also afect the time-step through its effect on the elastic forces?

In the following link you can find the results of a simple test showing this error and a video: results

Here you can see my repository with the branch containing the modifications: branch

Please, let me know what you think about this!

Cheers,
David

10

Hi Magali,

This is all very interesting! I have been also working with some visco-elastic models and it has been a bit tricky. Until now, in my models I have been setting the elastic time-step to match the numerical time-step (In the recent tests using 1000 yrs for both). Increasing the elastic time-step will probably help me reducing surface deflections. Have you found problems or issues using a smaller numerical time-step than the elastic time-step?

Also, the elastic time-step chosen will produce very significant changes in the results. The approach would be to tune the elastic time-step to produce velocities that match the time-scales of deformation that I expect? Or is there a rule of thumb to choose a “fair” elastic time-step?

Cheers,
David

11

Hello David,

I carefully read the Farrington et al Gcubed, 2014 paper and the original Moresi et al 2003 paper that describe the
visco-elastic method that is implemented in Aspect. Here’s a paragraph from Farrington about t_c the computational
time-step and t_e the elastic time-step:

A decoupling of the observation time and numerical time step is required [Moresi et al., 2003]. This can be achieved by introducing an elastic time step, te, that ensures the required observation time is resolved, distinct from the smaller numerical time step, tc . Since te > tc , we must store the stress history over several time steps…

Also, in one of the benchmark cases in the Farrington paper, she sets t_c = t_e/3. In this example, (torsion shear test), as t_e is reduce the peak elastic stress is recovered more accurately, but further reducing it, has little effect.

So, from reading these papers, I take it that it that t_e is required to be greater than t_c AND stress averaging should always be on.

From reading the Moresi paper and another reference in that paper (by Funicello), it appears that the elastic time-step should really be thought of the time over which you expect to see visco-elastic relaxation in the region of interested (e.g., the bending hinge of a slab). A first estimate of t_e is the Maxwell time t_m = eta/mu for the deforming layer (so in the subduction case, its the maxwell time for the most viscous part of the slab, but not the much weaker crust or the mantle).
I still need to test how sensitive the overall dynamics is to t_e.

I don’t know what kind of problem you working on, one in which you are interested in very short time-scale elastic response (like deformation related to faulting in the crust), which should have a quite short elastic time-scale, or more
the “long-term” tectonics time-scale (like the bending plate), which case a longer elastic time scale is appropriate.

I’m happy to discuss more, this is a little bit tricky and I’m still figuring this out also.
Magali

12

Hi Magali,

Thanks a lot for your comments, this has been really helpful! I will do more tests with this because I am interested in “long term” tectonics so I definitely need to increase a lot the elastic time-step.

Cheers,
David

13

Hi David,

Apologies for the delay in getting back to this discussion. At first glance, your code looks correct. Up to time step 19 (last output), are the viscosity values you obtained consistent with the analytical formulation?

Also, would it be possible for you to re-run this model and output every time step? The last output was at time step 19, but the models continues on for a bit, even though the time step is effectively approaching 0.

This can be achieved by setting:
set Time between text output = 0

Also, can you add material statistics to the list of post processors? This will output the effective viscosity to the log.txt and statistics file, which will make it a bit easier to plot this value against the analytical solution over time.

Thanks again sorry again for the delay in getting back to this discussion!

Cheers,
John

14

Hi John,

Thanks a lot for your reply. I re-ran the model with a lower viscosity to see the approach of the effective viscosity to the true viscosity before the drop in the time-step (You can see a plot of this in the folder). I also attached the results with output at every time-step and with material statistics.

What do you think about this?

Cheers,
David

15

Hi David,

Thanks for uploading the new solution results. However, I am still a bit confused about the source of the error.

What is happening is you are getting localized viscoelastic stress accumulation, which eventually leads to a propagating runaway effect. However, I am a bit confused as to why this is happening as there are no boundary driving forces or internal buoyancy variations. So, it looks like viscoelastic stresses initially develop from numerical noise and then grow due to the new formulation suggested here.

For the next step, I suggest trying your modified code with one of the viscoelastic benchmarks that has an analytical solution (viscoelastic_stress_build-up) or the bending beam problem. In particular, I think you should be able to derive a new analytical solution for the viscoelastic_stress_build-up problem based on your modified formulation.

Does this sound like a reasonable plan?

Cheers,
John

16

Hi John,

Thanks a lot for your reply and suggestion! This sounds like a good plan. I will use the viscoelastic_stress_build-up benchmark to do some tests with the code.

Cheers,
David