I am running some crustal models implementing visco-plastic rheology. During my runs, regardless of the adopted crustal geometry/density/viscosity, I have noticed that the total pressure is negative at shallow depths. I have incorporated a free surface at the top boundary, free slip at the bottom boundary and a very low, symmetric extensional velocity at the lateral boundaries. At first, I thought that this result comes from the dynamic pressure component caused by extension. However, the dynamic pressure, with velocities of a few mm/yr, should be very low, and should not exceed the lithostatic pressure, except maybe at the surface. Maybe you can explain to me why this happens?
Thank you very much!
Thank you for posting this question to the forum! Short answer, negative pressures near the model surface is quite common in extension models, even for slow extension rates. As you stated, this arises from the dynamic pressure component in extension.
How far down (depth) do the negative pressures extend in your models?
For a lithospheric-scale model with resolution on the order of 0.5-2 km, my recollection is that the negative pressures would normally be in the top 1 or 2 elements.
For reference, the default in the visco plastic model is to set the pressure for the Drucker Prager plasticity to 0 if it is negative:
Please let us know if you would like to discuss further, particularly if you think the model results are show unrealistically large regions with negative pressures.
Thank you very much for your quick reply. I have made sure that negative pressures are not allowed for the Drucker-Prager plasticity, but I am trying to understand why this problem arises. Right now, I have a negative pressure in the top 2 km of my model, corresponding to 4 cells with a dimension of 500m. Maybe my estimation is completely wrong, but by applying an extensional velocity of 1 mm/yr, being the dynamic pressure P_d = 0.5rhou^2, I should have values of P_d in the order of 10^-17 – 10^-18 Pa. If this is true, negative values of total pressure should never appear, except at the surface where P=0. However, when I check the output values from ASPECT for the nonadiabatic pressure, they range between 10^6 and 10^8 Pa. Can you give me some explanation about this difference in estimations?
@Alessio The fundamental problem is that we think that the pressure at the surface should be zero, but that’s not what the mathematics provide. The Stokes equations produce a pressure that is only determined up to an additive constant. So we can fix the pressure at a point, or we can select the average pressure along the surface to be zero, or we can select the average pressure throughout the domain to be zero. All of these are valid choices for the Stokes equations, and if you look at your input file (or the
parameters.prm file in your output directory), I suspect that you have the second of the options above selected. But that only means that the average surface pressure is zero, not that the pressure is zero along the entire surface.
Hi Alessio, Hi Wolfgang,
@Alessio to follow on @bangerth point above, we actually specifically disable pressure normalization with the free surface:
Regarding the values of the dynamic pressure, I think to first order the values you are seeing are not unreasonable.
If one just considers a strain rate of 1.e-16 (~ average value expected for 1 mm/yr over a 200 km wide model) and a viscosity of 1e22 or 1e24 Pa s, this gives stress value of 1 or 100 MPa (i.e, in the range you are seeing).
Is their localized deformation in your model (i.e., are there strong gradients in velocity and viscosity due to faulting)?
If yes, can you run the same model with a fixed viscosity for one time step and see if the predicted stresses and dynamic pressure is consistent with the strain rate and viscosity value?
Dear John and Wolfgang,
thank you very much for your fast and thorough reply. Actually, I have not tried to run a model with fixed viscosity, and it is indeed a test that I should carry out. I will do that and keep you posted.
Also, thank you very much for your suggestion to disable the pressure normalization with the free surface. However, in the prm file I already disabled the pressure normalization (set Pressure normalization = no)
I will keep on posting my results on this thread.