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

# 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