ASPECT: problem with the free surface


#1

Dear all,

I have a problem using the free surface boundary condition.
I have a box domain characterized by mantle and crust materials. In particular, the crustal level is defined by two materials, oceanic crust to the left and continental crust to the right, with different densities (2900 and 2800 Kg/m3, respectively) but same thickness. Left, right and bottom boundaries are free slip and the top boundary is free surface. No extra forces or velocities have been imposed. Using these conditions the model is not stable and it results in large vertical velocities in the mantle (3 m/yr, see attached picture), while I would expect velocities very small. I also copy here the input file (free_surface_ocean_continent.prm).

When I change the thickness of the two crusts (i.e. 10 km for the oceanic crust and 30 km for the continental crust) the vertical velocities dramatically increase up to 20 m/yr!

When I use the same density and thickness for the crust (basically a homogeneous crustal level), the model works fine, and the resulting velocities are very small (10^-9 m/yr).

It looks like that the free surface is not stable when different densities occur along the top boundary.

Do you have an idea about the problem with the free surface?

Many thanks in advance.
Regards,
Manuel

INPUT FREE SURFACE OCEAN-CONTINENT

Global parameters

set Dimension = 2
set Start time = 0
set End time = 5e6
set Use years in output instead of seconds = true
set Nonlinear solver scheme = single Advection, iterated Stokes
set Nonlinear solver tolerance = 1e-4
set Max nonlinear iterations = 10
set CFL number = 0.5
set Output directory = output-continental_extension
set Timing output frequency = 1
set Pressure normalization = no

Parameters describing the model

subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 100
set Y repetitions = 50
set X extent = 1400e3
set Y extent = 700e3
end
end

subsection Free surface
set Free surface boundary indicators = top
set Surface velocity projection = vertical
end

subsection Boundary velocity model
set Tangential velocity boundary indicators = bottom, left, right
end

subsection Compositional fields
set Number of fields = 4
set Names of fields = upper, lower, mantle, ocean
end

subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = if(x>=700e3 && y>=680.e3, 1, 0);
if(x>=700e3 && y<680e3 && y>=670.e3, 1, 0);
if(y<670.e3,1, 0);
if(x<700e3 && y>=670.e3, 1, 0);
end
end

subsection Boundary composition model
set Fixed composition boundary indicators = bottom
set List of model names = initial composition
end

subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box

subsection Box
set Bottom temperature = 1573
set Top temperature = 273
end
end

subsection Initial temperature model
set Model name = function
subsection Function
set Variable names = x,y
set Function constants = h=700e3,ts1=273,ts2=681.5714,ts3=823.,
k1=2.5,k2=2.5,k3=3.3,A=1.5e-6,
qs1=0.0653571,qs2=0.035357,qs3=0.035357,qb3=0.035357
set Function expression = if( (h-y)<=20.e3,
ts1 + (qs1/k1)(h-y) - (A(h-y)(h-y))/(2.0k1),
if((h-y)>20.e3 && (h-y)<=30.e3, ts2 + (qs2/k2)(h-y-20.e3),
if((h-y)>=30.e3 && (h-y)<100.e3, ts3 + (qs3/k3)
(h-y-30.e3),1573)));
end
end

subsection Heating model
set List of model names = adiabatic heating, shear heating
end

Material model

Rheology: Non-linear viscous flow and Drucker Prager Plasticity

subsection Material model
set Model name = visco plastic

subsection Visco Plastic

# Reference temperature and viscosity
set Reference temperature = 293
set Reference viscosity = 1e22

set Minimum strain rate = 1.e-20
set Reference strain rate = 1.e-16

set Minimum viscosity = 1e18
set Maximum viscosity = 1e26

set Thermal diffusivities = 1.737373e-6, 0.965517e-6, 0.965517e-6, 1.737373e-6, 0.858086e-6
set Heat capacities       =       1000.,        750.,        750.,       1000.,        750.
set Densities             =        3300,        2800,        2900,        3300,       3000.
set Thermal expansivities =        3e-5,        3e-5,        3e-5,        3e-5,        3e-5

set Viscosity averaging scheme = harmonic

set Viscous flow law = dislocation

set Prefactors for dislocation creep          = 2.41e-16, 1.13e-28, 1.13e-28, 2.41e-16, 1.99e-11
set Stress exponents for dislocation creep    =      3.5,      3.2,      3.2,      3.5,      3.4
set Activation energies for dislocation creep =   540.e3,   123.e3,   123.e3,   540.e3,   497.e3
set Activation volumes for dislocation creep  =   15.e-6,       0.,       0.,   15.e-6,       0.

# Plasticity parameters
set Angles of internal friction =   20.,   20.,  20.,    20.,   20.
set Cohesions                   = 20.e6, 20.e6, 20.e6, 20.e6, 20.e6

end
end

subsection Gravity model
set Model name = vertical

subsection Vertical
set Magnitude = 9.81
end
end

Post processing

subsection Postprocess
set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization
subsection Visualization
set List of output variables = density, viscosity, strain rate, error indicator
set Output format = vtu
set Time between graphical output = 1e6
set Interpolate output = true
end
end

subsection Solver parameters
subsection Stokes solver parameters
set Number of cheap Stokes solver steps = 0
end
end


#2

Hi Manuel,
great to hear that you are interested in using ASPECT.

About the problem you described, and just to see if I understand correctly: It seems you are aware that the imbalance in the crustal densities will lead to a velocity that will create some topography, however you assume the created velocities are bigger than what you expect? I agree that 3 m/s seems too large, but I have two questions about this:

  1. In the picture you attached I only see velocities of 3 cm/yr (3e-2), did you unintentionally post a different model?
  2. There could be an interaction between your visco plastic rheology and this density imbalance, i.e. it might not be the free surface, but the plastic failure that creates these large velocities. How does the viscosity distribution look like in your model, and what happens if you for example simplify the model to be isoviscous? I would guess this should resolve the instabilities and would lead us closer to the origin of these large velocities.

Hope that helps a bit. Let us know how it goes.

Best,
Rene


#3

Hi Manuel,

Echoing Rene, glad to hear you are trying ASPECT and thank you for posting on the forum.

I also agree with Rene that this behavior is not necessarily unexpected or unphysical given the initial density contrast and low viscosities that can result from the visco plastic rheology.

In oceanic-continental subduction models I often see this type of behavior (large vertical velocities) in the first few time steps when the model is isostatically adjusting. We dealt with this issues by either

  1. Using an isoviscous rheology until the model had isostatically adjusted (Rene’s suggestion) or

  2. Applying an initial topographic gradient to account for lateral density variations (e.g., enforce isostasy).

Cheers,
John


#4

Dear Manuel,

The reason for the fast velocities is because the model is not in isostatic equilibrium because the initial surface

height is flat. As you probably know from doing a simple isostasy calculation, the side of the model with the less

dense crust should have higher topography then the other side. When the model is in isostatic equilibrium the pressure below the layers with density anomalies will be the same. So, the first thing that happens in the simulation is there is flow drive by the different pressures below the crustal layer. However, this will lead to a problem known

as the “drunken sailor” phenomena in which the top surface will oscillate up and down. The paper by Kaus et al.,

attached describes this problem in more detail.

I am not familiar enough with the free-surface implementation in Aspect to answer as to what the solution to this problem - others on this e-mail list will have more information. One possible option is to start out by imposing an initial surface topography that is consistent with your density structure (calculate by hand what the isostatic topography should be). This should help to minimize the time it takes the simulation to get out of any transient behavior due to the initial set-up.

-Magali


#5

Many thanks for all replies!
I will try your suggestions and I let you know the result

Cheers
Manuel