# ASPECT: problem with the free surface

**URL:** <https://community.geodynamics.org/t/aspect-problem-with-the-free-surface/264>\
**Category:** ASPECT\
**Created:** [December 11, 2018, 3:30pm UTC](https://community.geodynamics.org/t/aspect-problem-with-the-free-surface/264 "2018-12-11T15:30:59Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![manuelroda](https://avatars.discourse-cdn.com/v4/letter/m/b4bc9f/32.png) [@manuelroda](https://community.geodynamics.org/u/manuelroda)\
**Post date:** [December 11, 2018, 3:30pm UTC](https://community.geodynamics.org/t/aspect-problem-with-the-free-surface/264/1 "2018-12-11T15:30:59Z")

</div>

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

 ![Vy_map](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/1X/9814e31333889248a709f3a8920fa59e52148cd6.png)

#### 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

---

<div class="post-metadata">

**Author:** ![gassmoeller](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/gassmoeller/32/11_2.png) [@gassmoeller](https://community.geodynamics.org/u/gassmoeller)\
**Post date:** [December 12, 2018, 2:05am UTC](https://community.geodynamics.org/t/aspect-problem-with-the-free-surface/264/2 "2018-12-12T02:05:33Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [December 12, 2018, 3:26am UTC](https://community.geodynamics.org/t/aspect-problem-with-the-free-surface/264/3 "2018-12-12T03:26:06Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mibillen](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/mibillen/32/169_2.png) [@mibillen](https://community.geodynamics.org/u/mibillen)\
**Post date:** [December 12, 2018, 8:11am UTC](https://community.geodynamics.org/t/aspect-problem-with-the-free-surface/264/4 "2018-12-12T08:11:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![manuelroda](https://avatars.discourse-cdn.com/v4/letter/m/b4bc9f/32.png) [@manuelroda](https://community.geodynamics.org/u/manuelroda)\
**Post date:** [December 12, 2018, 9:10am UTC](https://community.geodynamics.org/t/aspect-problem-with-the-free-surface/264/5 "2018-12-12T09:10:42Z")

</div>

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

Cheers  
Manuel