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 , 10cm/year, -1cm/year) ; 0;
!# set Function expression = if( z>=litho_d,5cm/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