Principal stress orientations in viscous vs visco-elastic models

Hi all,
We (@mibillen) have been looking at the principal stress visualization outputs in our 2D visco-elastic-plastic subduction models and were surprised by the principal stress orientations calculated in the slab. We ran a 2D visco-plastic subduction model for comparison, and there was significant (90 degrees in some areas) differences in principal stress orientations between a visco-elastic-model and a viscous model (with the viscous model showing the physically more expected stress orientations in the slab). The visco-elastic orientations seemed wrong, and we’re surprised, but if you’re not comparing viscous to visco-elastic models you might not notice. I have now run some small test models to illustrate what we were seeing and what seems like it might be incorrect.

I modified the bending beam box model to suspend the beam as a dense block of rock from the top of box. We would expect this beam to be in down-dip tension as gravity is pulling it downward. I ran two versions of this model that are identical except for the “Enable elasticity” formulation being set to true in one model, and false in the other. ( VEtopbeam.prm (6.6 KB), Vtopbeam.prm (6.6 KB))

These outputs show the results of running the two versions of the beam model. You can see on the left, in the viscoelastic model, that the sigma_1 direction indicates the beam is in down-dip compression, contrary to what we physically would expect. On the right, in the viscous model, the sigma_1 direction indicates the beam is in down-dip tension, as we would physically expect.

In aspect/source/postproccess/visualization/, we saw the note on compressive stress being positive followed by the flag for incorporating elasticity:

I tested changing all of the “+=” signs within this elastic part to “-=”, recompiled, and then re-ran the above models. The new “principal_stress_direction_1” orientations are rotated and now show down-dip tension in the hanging beam both in the viscoelastic and in the viscous model, as expected:

@jbnaliboff @anne-glerum We think this is the correct orientation and that this sign in incorporating the elasticity was just incorrect, unless there is something physical I am misunderstanding between the viscous and viscoelastic behavior?

Additionally, this is just in the post processing, but I was not sure where to find within the solvers how the stress tensor was being handled sign wise while solving the stokes system?

Maybe related?: fixed stress visualization for elastic materials by bobmyhill · Pull Request #4375 · geodynamics/aspect · GitHub

Hi Becky, thanks for reporting this bug. You are right, this looks incorrect and is likely caused by using the wrong stress in the postprocessor if elasticity is enabled. I had a look at your fix, but it doesnt look correct to me either. In your fix you first compute the full non-elastic stress, then subtract the full stress including elastic stresses, which should give you only the elastic stresses. Could you check if the fix I proposed here works for you: Fix elastic stress used in principal_stress postprocessor by gassmoeller · Pull Request #5480 · geodynamics/aspect · GitHub? This should use either the full stress including elasticity (if elasticity is on), or just the usual stress computed from the strain rate (if elasticity is off). I havent tested this, but it is essentially the same fix that Bob linked for another postprocessor.

Let me know, It would also be great if we could get the elastic beam example you used as a test case for this fix (currently this postprocessor option was untested, this was part of the problem of course).


Thank you Rene, We had come to the same conclusion after reading Bob’s pull request, but Becky had a job interview yesterday so was focused on preparing for that.

Becky will test this later this week, and will send you the test case so this can also be added to the tester.


Hi Rene,
Yeah as Magali replied, we see now that difference and looked at what Bob had linked. I tested the new PR and found that it alters the sigma1 direction in some areas, but that the viscoelastic stress orientations in the hanging beam are still in the opposite direction:

I think this has to do with the sign convention still. If I am understanding correctly, inside of the code, the sign convention follows the engineering principles of compression having a negative sign, and then in postprocessing, values are converted to the geologic convention of compression having a positive sign. In line 144, the stress sign is being switched for the viscous case:
Screen Shot 2023-11-09 at 10.11.02 AM
In the viscoelastic case though, the stress components are just being brought in as their existing values stored in the compositional fields, so it is still in the engineering sign convention from inside the code and the sign change during postprocessing is not being accounted for in these lines:

I made some commented suggestions to the code change in PR 5480 that I think would be correct now if I am understanding everything. I have tested this and it does produce downdip tension in the hanging beam test cases I have.

Also yes, I provided my two hanging beam prms (both the viscous and viscoelastic cases) above and am happy to contribute these as test cases!


Great catch Becky, thanks!

I’m having some problems with building ASPECT on my machine (xcode/mac problems that I’m too stupid to figure out). With that and usual lecturer busywork (which I’m also too stupid to figure out) I probably won’t get round to fixing PRs until December. You’re welcome to take over that PR if you want.