Asking for advice for boundary condtion for subduction problems

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

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

1 Like

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