# Asking for advice for boundary condtion for subduction problems

**URL:** <https://community.geodynamics.org/t/asking-for-advice-for-boundary-condtion-for-subduction-problems/658>\
**Category:** ASPECT\
**Created:** [June 9, 2019, 4:17am UTC](https://community.geodynamics.org/t/asking-for-advice-for-boundary-condtion-for-subduction-problems/658 "2019-06-09T04:17:09Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![sungho](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/sungho/32/365_2.png) [@sungho](https://community.geodynamics.org/u/sungho)\
**Post date:** [June 9, 2019, 4:17am UTC](https://community.geodynamics.org/t/asking-for-advice-for-boundary-condtion-for-subduction-problems/658/1 "2019-06-09T04:17:09Z")

</div>

Hello all.

Although my background is hydrogeology and geomechanics, I am a Ph.D. student studying recently geodynamics and seismology because I have a private circumstance.

Recently, I am trying to build a model for plate subduction model, however, I don’t understand exactly boundary conditions used in that kind model because my understanding of the conceptual model for global tectonics and geodynamics are poor.

I’ve read several papers related to geodynamics numerical modeling.  
The papers are saying that prescribe velocity condition is assigned on an oceanic plate (top boundary) or on oceanic crust + oceanic lithosphere (side boundary).

I think either is a “make sense” in the conceptual aspect.

However, I can’t figure out what is the difference between them in terms of numerical and physical things.

In my knowledge, they must show a similar phenomenon based on the assumption of the rigid plate, however, the lithosphere beneath oceanic crust assigned prescribed velocity in my model is shortened even though I give constant composition boundary condition on the side walls.

This is not what I want.

So, my point is asking for advice for some general conditions (boundary condition, non-linear property model ) for a subduction (slab) model.

I know this is an ambiguous question but I need a help.

Thanks

Sungho

 ![visit0000](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/1X/86d5efaca631cf1e5245f076722b5a3101707b1b.png)

# 2D, Dimensional subdcution model

set World builder file = flat5.wb  
set Output directory = output-flat-bc8

set Dimension = 2  
set Use years in output instead of seconds = true  
set Resume computation = false  
set Start time = 0  
set End time = 1e6  
set Use years in output instead of seconds = true  
set Nonlinear solver scheme = single Advection, iterated Stokes  
set Nonlinear solver tolerance = 1e-2  
set Max nonlinear iterations = 10  
set CFL number = 0.5  
set Timing output frequency = 1

set Pressure normalization = surface  
set Adiabatic surface temperature = 1600 # K

subsection Geometry model  
set Model name = box with lithosphere boundary indicators  
subsection Box with lithosphere boundary indicators  
set X extent = 2000e3  
set Y extent = 660e3  
set Lithospheric thickness = 100e3  
set X repetitions = 4  
end  
end

!# Advecting the free surface vertically rather than  
!# in the surface normal direction can result in a  
!# more stable mesh when the deformation is large  
!#subsection Free surface  
!# set Free surface boundary indicators = top  
!# set Surface velocity projection = vertical  
!#end

subsection Gravity model  
set Model name = vertical  
subsection Vertical  
set Magnitude = 9.81  
end  
end

subsection Compositional fields  
set Number of fields = 4  
set Names of fields = MT, OC, UC, LC  
end

subsection Initial composition model  
set Model name = world builder  
end

subsection Initial temperature model  
set Model name = world builder  
end

!# 0 left, 1 right, 2 bottom, 3 top, 4 left litho, 5 right litho  
!# top and bottom is not fixed, it could be changed next model  
subsection Boundary temperature model  
set Fixed temperature boundary indicators = 2, 3  
set List of model names = initial temperature  
end

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

subsection Boundary composition model  
set Fixed composition boundary indicators = 0, 1, 4, 5  
set List of model names = initial composition  
end

!# 0 left, 1 right, 2 bottom, 3 top, 4 left litho, 5 right litho

subsection Boundary velocity model  
set Tangential velocity boundary indicators = 0, 1, 4, 5  
set Prescribed velocity boundary indicators = 3: function  
set Zero velocity boundary indicators = 2  
subsection Function  
set Variable names= x,z  
set Function constants = cm=0.01, year=1, lith\_th=100e3, model\_height=660e3, litho\_d=560e3  
set Function expression = if (x\<=1000e3 , 10_cm/year, -1_cm/year) ; 0;  
!# set Function expression = if( z\>=litho\_d,5_cm/year_(z-litho\_d)/lith\_th, 0); 0  
end  
end

!### lithostatic pressure dose not be together free surface boundary

!#subsection Boundary traction model  
!# set Prescribed traction boundary indicators = bottom: initial lithostatic pressure #, 1: lithostatic pressure  
!#end

!#subsection Boundary traction model  
!# subsection Initial lithostatic pressure  
!# set Representative point = 0,0  
!# end  
!#end

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, 1.737373e-6, 0.858086e-6, 0.965517e-6, 0.965517e-6  
set Heat capacities = 1000., 1000., 750., 750., 750.  
set Densities = 3300, 3300, 3000., 2800, 2900  
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, 2.41e-16, 1.99e-11, 1.13e-28, 1.13e-28  
set Stress exponents for dislocation creep = 3.5, 3.5, 3.4, 3.2, 3.2  
set Activation energies for dislocation creep = 540.e3, 540.e3, 497.e3, 123.e3, 123.e3  
set Activation volumes for dislocation creep = 15.e-6, 15.e-6, 0., 0., 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 Mesh refinement  
set Initial global refinement = 3  
set Initial adaptive refinement = 4  
set Refinement fraction = 0.95  
set Coarsening fraction = 0.05  
set Time steps between mesh refinement = 10  
set Strategy= composition, temperature  
end

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 = 1e3  
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:** ![jdannberg](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/jdannberg/32/26_2.png) [@jdannberg](https://community.geodynamics.org/u/jdannberg)\
**Post date:** [June 17, 2019, 3:22pm UTC](https://community.geodynamics.org/t/asking-for-advice-for-boundary-condtion-for-subduction-problems/658/2 "2019-06-17T15:22:17Z")

</div>

Hi Sungho,

if I understand you correctly, you want the temperature at the left model boundary to be similar to the lithospheric temperature profile you start out with.  
In your input file, you only prescribe the temperature at the top and the bottom of your model domain, but you have inflow from the left. So you have to tell ASPECT what temperature it should assign to that inflowing material. You can do that by adding the corresponding boundary indicator to your ‘Fixed temperature boundary indicators’.

In your specific case, I believe that would be indicator 4 (for the lithosphere on the left model boundary).  
So

> set Fixed temperature boundary indicators = 2, 3, 4

Does that fix your problem?

Best regards,  
Juliane

---

<div class="post-metadata">

**Author:** ![sungho](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/sungho/32/365_2.png) [@sungho](https://community.geodynamics.org/u/sungho)\
**Post date:** [June 18, 2019, 7:41am UTC](https://community.geodynamics.org/t/asking-for-advice-for-boundary-condtion-for-subduction-problems/658/3 "2019-06-18T07:41:46Z")

</div>

Hi Juliane

Thank you for your comment.  
However, not really is it work.

Fixed temperature on the left boundary stay only temperature on the thin left vertical layer.  
Right after the left boundary, shortening is still overwhelming.

It would seem related to prescribed velocity on the surface is too high compared to heat conduction rate.

But, I’m not sure 10 cm/year is high velocity or not comparing to heat conduction in this model.

Best regards  
Sungho