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