Hi everyone,
I am working with the viscoplastic material model and am trying to simulate plate-like behavior in the presence of highly viscous continental material. My objective is to assign a much higher viscosity to the continents than to the surrounding mantle. However, I have noticed that regardless of how I modify the continental rheology parameters (e.g., prefactors, activation energies, or activation volumes), the effective viscosity of the continents never exceeds roughly twice the mantle viscosity. This behavior appears to be independent of the specific parameter values used.
Interestingly, this issue only arises when the continental viscosity is prescribed to be higher than the mantle viscosity. If I instead make the continents weaker than the mantle, I can obtain essentially any viscosity contrast without difficulty.
The behavior is particularly clear when using a P–T-independent viscoplastic rheology (i.e., setting E=0 and V=0), where a constant viscosity contrast can be directly reflected in the solution. I have attached the relevant .prm file for reference.
Could this be related to viscosity averaging, compositional field averaging, or some other aspect of the viscoplastic material model? More generally, is it possible within the viscoplastic framework to make continental material sufficiently strong that its viscosity effectively reaches the global upper viscosity cutoff?
One possible workaround would be to prescribe the continents through a World Builder file, but I would prefer to avoid that approach since I need the continental structure to persist naturally through model restarts.
I would greatly appreciate any suggestions or insights. I am happy to provide any additional information.
Thank you very much for your help.
Best,
Arnab
aspect_model_vp.prm (8.7 KB)
Hi Arnab,
Thanks for posting on the forum! This does sound like it could be an issue with the averaging, what is the resolution of your mesh and what is the thickness of your continents? It would be helpful to see your continent composition with your mesh overlain in Paraview. I suspect that your cell size is roughly equal to the thickness of your continents, which combined with the harmonic viscosity averaging scheme you are using would explain why you can reduce the viscosity contrast, but not increase it.
Harmonic averaging sums up the viscosity of the different compositions in parallel, which will always bias the result to the lower viscosity value. So in the case of a continent viscosity (eta_c) and a mantle viscosity (eta_m), your harmonically averaged viscosity would be:
1/eta = (1/eta_m + 1/eta_c)
eta = eta_m * eta_c / (eta_m + eta_c)
Playing with the values of eta_m and eta_c, you can convince yourself that the viscosity will tend towards the lower of the two. Looking forward to hearing if this helps address your issue.
Cheers,
Daniel
Hi Arnab,
I agree with Daniel. Averaging might be the issue here. Looking at the input file, it seems that, at the location where you define your continental (your 2nd composition) to be 1, i.e. ‘true’, you have also defined mantle material (your 1st composition) to be 1, because mantle material is set to 1 everywhere:
set Function expression = 1; if(r>6171000 && phi>=0 && phi<1.5707, 1, 0)
So the viscosity at those locations will be harmonic average (as you set that in the input file), which, as Daniel describes will tend to be closer to the weaker of the two.
I wonder if you need ‘mantle’ material at all in your models. You could just define continental material as the only composition, and wherever you don’t have continental material, the ‘default’ value will be taken, which you defined as mantle material already. Alternatively, set mantle material to be 1 everywhere EXCEPT where the continent is located.
Cheers,
Jeroen