Issues with pressure on top surface

ASPECT
1

Hi all,

I have always had an issue with pressure on the top surface.

How can I make the pressure at a free upper surface to have a specified value?

I have tried:

a) Setting normalization

set Pressure normalization = surface
set Surface pressure = 1e5

for example if I want a typical atmospheric pressure on the top surface of the model. I guess in this case Aspect integrates the Stokes eqs., and finally at the moment of postprocessing it adds a constant so that the pressure field starts at my desired value.

b) Prescribing a traction BC at that upper surface

subsection Boundary traction model
set Prescribed traction boundary indicators = top: function

subsection Function
set Coordinate system = cartesian
set Variable names = x,y,z
set Function constants = Ptop=1e5
set Function expression = 0 ; 0 ; Ptop

 ##### I also tried Function expression = 0 ; 0 ; -Ptop    to account for the vector nature

end

end

I guess in these cases Aspect integrates the component term dp/dz using that prescribed value equal to ±Ptop.

c) I also tried setting (a) and (b) simultaneously.

But when I plot the model results in Paraview, I don’t get the desired value of positive pressure 1e5 Pa. Usually I get a negative pressure (which is not exactly equal to -Ptop).
I totally know this is not a problem, Aspect does calculate the pressure field inside correctly, as on the plots it varies as it should. Simple calculus tells that when integrating the Stokes flow equation, the pressure needs a constant if we want absolute values.

But sometimes, I would like to have a pressure plot, and having negative values on top looks odd and reckless.
Most importantly, sometimes I certainly want my model to actually have a pressure on top ( I think in some cases, if internal stresses are not too high, and if the model runs for a long time, that surface pressure may potentially be a small controlling factor).

I have tried this in several models, sometimes also with a loaded topography. The results are always similar: pressure field has the correct variation, it was seemingly well calculated but without the integration constant as what I see on top is never the value I want.

It may be an error in Paraview when reading the data…

What should I do?

cheersmanycheers,
Felipe

2

Hi Felipe,

There are multiple ways to prescribe a traction on the top boundary. You noted one method with the set Pressure normalization and set Surface Pressure parameters.

The value you selected (1e5) is only 0.1 MPa, and thus negative non-lithostatic pressures can easily arise. You could use a value of 1 or 10 MPa (depending on magnitude of negative surface pressures), which is equivalent to ~ 100 or 1000 m of water overburden.

Using a boundary traction on the upper surface also allows you to prescribe a surface overburden, but of course changes the dynamics of deformation near the surface as material can pass through that boundary.

There are also examples of how to combine a free surface with a surface boundary traction, which match analytical solutions:

So, my suggestion is to follow one of these three methods and just prescribe a higher surface traction (pressure) if you don’t want to see negative pressures in your plots. Start with the examples linked above, as those have been well documented.

John