Hi everyone,
I am working on a continental collision model where an oceanic plate subducts beneath a continent while dragging another continental piece which eventually results in a continental collision. The model seems to be running fine but I am getting this grid-like pattern in the viscosity field after the initial timestep which makes me believe something is wrong. See the figure for clarification.
Example viscosity field:
I am using Geodynamic World Builder to create geometry. 6 km thick oceanic crust on the subducting plate is modified to act as a lubricant layer (I decreased the viscosity of this layer artificially to create sufficient lubrication). You can find .prm and .wb files below. Other than the examples below; I have tried to switch off the free surface, I have tried models with just diffusion creep and composite flow laws but the results are similar to the example screenshots.
.prm file:
#### World Builder parameters which can be used by ASPECT ######## World Builder parameters which can be used by ASPECT ####
set World builder file = sub2col_base_2PlatDip45.wb
set Dimension = 2
set Start time = 0
set End time = 50e6
set Pressure normalization = no
set Maximum time step = 1e4
set Maximum first time step = 100
set Maximum relative increase in time step = 20
set Use years in output instead of seconds = true
set Nonlinear solver scheme = iterated Advection and Stokes
set Max nonlinear iterations = 30
set Nonlinear solver tolerance = 1e-3
set CFL number = 0.1
subsection Solver parameters
subsection Stokes solver parameters
set Linear solver tolerance = 1e-6
set Number of cheap Stokes solver steps = 100
end
end
set Output directory = sub2col_base_v31
subsection Initial temperature model
set Model name = world builder
end
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 1600
set Top temperature = 273
end
end
subsection Initial composition model
set Model name = world builder
end
#### parameters needed by ASPECT when using compositions ####
subsection Compositional fields
set Number of fields = 5
set Names of fields = Oceanic_Plate, Mantle, Plateau_Crust, Plateau_Mantle, Oceanic_Crust
end
subsection Geometry model
set Model name = box
subsection Box
set X extent = 2640e3
set Y extent = 660e3
# set Z extent = 660e3
set X repetitions = 8
set Y repetitions = 2
# set X periodic = true
end
end
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 5
end
##Free Surface
subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface
subsection Free surface
set Free surface stabilization theta = 0.6
set Surface velocity projection = vertical
end
end
##Boundary classifications
subsection Boundary velocity model
set Tangential velocity boundary indicators = bottom, left, right
end
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
set Viscous flow law = dislocation
# Reference temperature and viscosity
set Reference temperature = 293
set Reference viscosity = 1e22
# The reference strain-rate is used on the first non-linear iteration
# of the first time step when the velocity has not been determined yet.
set Reference strain rate = 1.e-16
set Minimum strain rate = 1.e-20
# Limit the viscosity with minimum and maximum values
set Minimum viscosity = 1e18
set Maximum viscosity = 1e23
# Material properties
set Thermal diffusivities = 1.333333e-6
set Heat capacities = 750.
#Oceanic_Plate, Mantle, Plateau_Crust, Plateau_Mantle, Oceanic_Crust
set Densities = 3300, 3300, 3300, 2900, 3300, 3300
set Thermal expansivities = 2e-5
# Harmonic viscosity averaging
set Viscosity averaging scheme = harmonic
set Prefactors for dislocation creep = 6.52e-16, 6.52e-16, 6.52e-16, 8.57e-28, 6.52e-16, 5.5e-18,
set Stress exponents for dislocation creep = 3.5, 3.5, 3.5, 4.0, 3.5, 3
set Activation energies for dislocation creep = 530.e3, 530.e3, 530.e3, 223.e3, 530.e3, 345.e3
set Activation volumes for dislocation creep = 18.e-6, 18.e-6, 18.e-6, 0., 18.e-6, 38.e-6
# Plasticity parameters
set Angles of internal friction = 5.7392, 5.7392, 5.7392, 5.7392, 5.7392, 5.7392
set Cohesions = 1e6
set Maximum yield stress = 4e8
end
end
subsection Gravity model
set Model name = vertical
end
# Post processing
subsection Postprocess
set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization, topography, dynamic topography
subsection Visualization
set List of output variables = density, viscosity, strain rate, error indicator
set Output format = vtu
set Time between graphical output = 5e5
end
subsection Topography
set Output to file = true
set Time between text output = 5e5
end
end
.wb file:
{
"version":"0.4",
"cross section":[[0,330e3],[2640e3,330e3]],
"coordinate system":{"model":"cartesian"},
"features":
[
{"model":"oceanic plate", "name":"Subducting Plate", "max depth":100e3,
"coordinates":[[0,0],[0,660e3],[1660e3,660e3],[1660e3,0e3]],
"temperature models":
[
{"model":"plate model", "max depth":100e3, "spreading velocity":0.01,
"ridge coordinates":[[0e3,-1],[0,2000e3]]}
],
"composition models":[
{"model":"uniform", "compositions":[4], "min depth":0, "max depth":6e3},
{"model":"uniform", "compositions":[0], "min depth":6e3,"max depth":100e3}]
},
{"model":"continental plate", "name":"Plateau_2", "coordinates":[[800e3,0],[800e3,660e3],[1200e3,660e3],[1200e3,0]],
"temperature models":[{"model":"linear", "max depth":120e3, "bottom temperature":1600}],
"composition models":[
{"model":"uniform", "compositions":[2], "min depth":0, "max depth":20e3},
{"model":"uniform", "compositions":[3], "min depth":20e3,"max depth":120e3}]
},
{"model":"continental plate", "name":"Plateau_1", "coordinates":[[2640e3,0],[2640e3,660e3],[1660e3,660e3],[1660e3,0]],
"temperature models":[{"model":"linear", "max depth":120e3, "bottom temperature":1600}],
"composition models":[
{"model":"uniform", "compositions":[2], "min depth":0, "max depth":20e3},
{"model":"uniform", "compositions":[3], "min depth":20e3,"max depth":120e3}]
},
{"model":"mantle layer", "name":"upper mantle", "min depth":100e3, "max depth":660e3,
"coordinates":[[0,0],[0,660e3],[2640e3,660e3],[2640e3,0]],
"temperature models":[{"model":"uniform", "temperature":1600}],
"composition models":[{"model":"uniform", "compositions":[1]}]
},
{"model":"subducting plate", "name":"Slab",
"coordinates":[[1660e3,660e3],[1660e3,0]],
"dip point":[1e7,-1e7],
"segments":[
{"length":250e3, "thickness":[100e3], "angle":[0,30]},
{"length":250e3, "thickness":[100e3], "angle":[30,45]}],
"temperature models":[{"model":"linear", "max distance slab top":100e3, "top temperature":273, "bottom temperature":1600}],
"composition models":[
{"model":"uniform", "compositions":[4], "min distance slab top":0, "max distance slab top":6e3},
{"model":"uniform", "compositions":[0], "min distance slab top":6e3,"max distance slab top":100e3}] }
]
}
