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.0*k1),

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