Negative pressures at surface

Alessio · May 19, 2023, 2:14pm

Hi all,
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!
Alessio

jbnaliboff · May 19, 2023, 5:54pm

Hi Alessio,

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:

github.com

geodynamics/aspect/blob/main/source/material_model/rheology/visco_plastic.cc#L291


      
                                                                    in.position[i]);
          output_parameters.current_friction_angles[j] = current_friction;
          output_parameters.current_cohesions[j] = current_cohesion;
          
          
// Step 5: plastic yielding
          
          
// Determine if the pressure used in Drucker Prager plasticity will be capped at 0 (default).
          // This may be necessary in models without gravity and when the dynamic stresses are much higher
          // than the lithostatic pressure.
          
          
double pressure_for_plasticity = in.pressure[i];
          if (allow_negative_pressures_in_plasticity == false)
            pressure_for_plasticity = std::max(in.pressure[i],0.0);
          
          
// Step 5a: calculate Drucker-Prager yield stress
          const double yield_stress = drucker_prager_plasticity.compute_yield_stress(current_cohesion,
                                                                                     current_friction,
                                                                                     pressure_for_plasticity,
                                                                                     drucker_prager_parameters.max_yield_stress);
          
          
// Step 5b: select if yield viscosity is based on Drucker Prager or stress limiter rheology

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.

Cheers,
John

Alessio · May 20, 2023, 10:20am

Dear John,
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?
Thanks again,
Alessio

bangerth · May 21, 2023, 10:25pm

@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.

Best
W.

jbnaliboff · May 22, 2023, 4:38pm

Hi Alessio, Hi Wolfgang,

@Alessio to follow on @bangerth point above, we actually specifically disable pressure normalization with the free surface:

github.com

geodynamics/aspect/blob/main/source/mesh_deformation/free_surface.cc#L41


      
          
          

          
namespace aspect
          {
            namespace MeshDeformation
            {
              template <int dim>
              void
              FreeSurface<dim>::initialize ()
              {
                // Pressure normalization doesn't really make sense with a free surface, and if we do
                // use it, we can run into problems with geometry_model->depth().
                AssertThrow ( this->get_parameters().pressure_normalization == "no",
                              ExcMessage("The free surface scheme can only be used with no pressure normalization") );
          
          
      // Check that we do not use the free surface on a boundary that has zero slip,
                // free slip or prescribed velocity boundary conditions on it.
          
          
      // Get the zero velocity boundary indicators
                std::set<types::boundary_id> velocity_boundary_indicators = this->get_boundary_velocity_manager().get_zero_boundary_velocity_indicators();

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?

Cheers,
John

Alessio · May 23, 2023, 8:37am

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.
Very best,
Alessio

Topic		Replies	Views
"Shrinking earth" from dynamic topography plugin ASPECT	33	915	April 12, 2023
Pressure on the surface for case with free surface boundary condition ASPECT	0	218	October 29, 2021
Left and Right Boundary producing different stress values for Continental extension ASPECT	2	135	November 11, 2023
A question about stress change due to inward velocity and outward velocity ASPECT	1	281	October 17, 2019
Problem with running 3D continental extension models ASPECT	10	187	March 28, 2024

Negative pressures at surface

Related topics