Hi all,
I have a subduction model I have been working on and I am trying to recreate it by using World Builder for the ease of future development, but I am running into this issue where I can’t get the same viscosity values even though I am using the same values for prefactor, activation energy, and activation volume etc. See the figure below for viscosity fields and viscosity change with the depth of the two models at x = 2000km, and for the prm files and the wb file. The bottom model is the one that utilizes World Builder.
Normal subduction model prm:
#### Global parameters
set Dimension = 2
set Start time = 0
set End time = 50e6
set Use years in output instead of seconds = true
set Nonlinear solver scheme = iterated Advection and Newton Stokes
set Nonlinear solver tolerance = 1e-3
set Maximum time step = 1e3
set Maximum first time step = 100
set Maximum relative increase in time step = 20
set Max nonlinear iterations = 30
set CFL number = 0.1
set Output directory = sub2col_v99
set Pressure normalization = no
subsection Solver parameters
subsection Stokes solver parameters
set Linear solver tolerance = 1e-4
set Number of cheap Stokes solver steps = 100
end
subsection Newton solver parameters
set Max pre-Newton nonlinear iterations = 100
set Nonlinear Newton solver switch tolerance = 1e-4
set Max Newton line search iterations = 0
set Maximum linear Stokes solver tolerance = 1e-1
set Use Newton residual scaling method = false
set Use Newton failsafe = true
set Stabilization preconditioner = SPD
set Stabilization velocity block = SPD
set Use Eisenstat Walker method for Picard iterations = true
end
end
# Model geometry (660 km deep aspect ratio 7:1)
subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 14
set Y repetitions = 2
set X extent = 4620e3
set Y extent = 660e3
set X periodic = true
end
end
# Mesh refinement specifications (no mesh refinement)
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 5
end
# Number and names of compositional fields
# crust = oceanic crust = weak layer on subducting plate
# plateau = buoyant continental crust that forms continental
subsection Compositional fields
set Number of fields = 2
set Names of fields = crust, plateau
end
# Spatial domain of different compositional fields
# For symbols used, see temperature distribution.
# Composition distribution pseudo code:
# For crust distribution:
# if (z<dc) then
# top layer: assign crust everywhere
# else
# if (distance from point (xt,zt) smaller than r and larger than r-dc,
# and x>xt and z<zs ) then
# assign crust
# else
# no crust
# endif
# endif
# For plateau distribution:
# if (x<xt, i.e. in subducting plate) then
# if (xpr-wp < x < xpr and z<dp) then
# assign plateau
# else
# no plateau
# endif
# else
# if (distance from point (xt,zt) larger than r and
# x<xt+wp and z<dp) then
# assign plateau
# else
# no plateau
# endif
# endif
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function constants = h=660e3, ar=7, ts=273, tm=1300, dc=6e3, xt=2640e3, wocrhs=1320e3, \
xpr=1000e3, wp=400e3, dp=40e3, \
r=6e5, zs=3e5, kappa=1e-6, aget=2.52e15
set Function expression = \
if( (h-y)<dc, \
1, \
if( sqrt((x-xt)*(x-xt)+(h-r-y)*(h-r-y))<r && \
sqrt((x-xt)*(x-xt)+(h-r-y)*(h-r-y))>(r-dc) && \
x>xt && \
y>h-zs, 1, 0) \
); \
if( x<xt, \
if( (x>xpr-wp && x<xpr && (h-y)<dp), 1, 0), \
if( sqrt((x-xt)*(x-xt)+(h-r-y)*(h-r-y))>r && \
(x<(xt+wp) && (h-y)<dp), 1, 0) \
)
end
end
# Temperature boundary conditions
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 1573
set Top temperature = 273
end
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 = normal
end
end
# Boundary classifications
subsection Boundary velocity model
set Tangential velocity boundary indicators = bottom
end
subsection Nullspace removal
set Remove nullspace = net x translation
end
# Initial temperature field
subsection Initial temperature model
set Model name = function
subsection Function
set Variable names = x,y
# h = box height
# ar = aspect ratio of box
# ts = surface temperature
# tm = mantle temperature
# dc = crustal thickness (= thickness of weak zone)
# xt = x-position of trench
# wocrhs = width of oceanic lithoshere on rhs of overriding continent
# xpr = x-coord of right edge of buoyant plateau
# wp = width of buoyant plateau
# dp = crustal thickness in buoyant plateau
# eps = a very small number to avoid division by zero age
# r = slab bending radius
# zs = max initial slab depth
# kappa = thermal diffusivity
# aget = age of plate, e.g. 100 (Myrs) *1e6*365*24*3600=100*3.15e13=3.15e15 sec
#
# Temperature distribution pseudo code:
# if (x<xt) then
# apply halfspace cooling (HSC) geotherm with age increasing from t=0 at lhs
# to age = aget at x=xt (the trencha location)
# else
# if (distance from slab curvature origin (at coordinates (x,z)=(xtr,r)) < r AND z<zs) then
# location is inside bending slab, so apply HSC geotherm with T=ts at radius r, and
# T increasing for radius smaller than r.
# else
# apply overriding plate linear geotherm:
# if x<ar*h-wocrhs (i.e. distance from rhs boundary > wocrhs) then
# apply HSC geotherm for plate age aget
# else
# apply HSC geotherm with age decreasing from
# age = aget at distance wocrhs from rhs boundary
# to age = 0 at rhs boundary
# endif
# endif
# endif
# with half-space geotherm: T=Tm*(1-erfc(z/(2*sqrt(kappa*t)))).
set Function constants = h=660e3, ar=7, ts=273, tm=1300, dc=6e3, xt=2640e3, wocrhs=1320e3, \
xpr=1000e3, wp=400e3, dp=40e3, eps=1e-5,\
r=6e5, zs=3e5, kappa=1e-6, aget=2.52e15
set Function expression = \
if(x<xt, \
ts+tm*(1-erfc((h-y)/(abs(2*sqrt(kappa*x/xt*aget)+eps)))), \
if((sqrt((x-xt)*(x-xt)+(h-r-y)*(h-r-y))<=r && y>h-zs), \
ts+tm*(1-erfc((r-sqrt((x-xt)*(x-xt)+(h-r-y)*(h-r-y)))/(2*sqrt(kappa*aget)))),\
if (x<ar*h-wocrhs, \
ts+tm*(1-erfc((h-y)/(2*sqrt(kappa*aget)))), \
ts+tm*(1-erfc((h-y)/(abs(2*sqrt(kappa*(ar*h-x)/wocrhs*aget))+eps))) \
) \
) \
)
end
end
# Constant internal heat production values (W/m^3)
subsection Heating model
set List of model names = compositional heating
subsection Compositional heating
set Compositional heating values = 0.
end
end
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
set Viscous flow law = composite
# 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
# 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.
set Densities = 3300, 3300, 2900
set Thermal expansivities = 2e-5
# Harmonic viscosity averaging
set Viscosity averaging scheme = harmonic
set Prefactors for diffusion creep = 1.5e-16, 1e-50, 1.5e-16
set Activation energies for diffusion creep = 375.e3, 0.e3, 375.e3
set Activation volumes for diffusion creep = 3.e-6, 0., 3.e-6
set Prefactors for dislocation creep = 1.1e-16, 1e-50, 1.1e-16
set Stress exponents for dislocation creep = 3.5, 1.0, 3.5
set Activation energies for dislocation creep = 530.e3, 0.0, 530.e3
set Activation volumes for dislocation creep = 14.e-6, 0.0, 14.e-6
# Plasticity parameters
set Angles of internal friction = 5.7392, 0., 5.7392
set Cohesions = 1e6, 5e6, 1e6
set Maximum yield stress = 4e8
end
end
# Gravity model
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 9.81
end
end
# Post processing
subsection Postprocess
set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization, topography
subsection Visualization
set List of output variables = density, viscosity, strain rate, error indicator, principal stress, shear stress
set Output format = vtu
set Time between graphical output = 1e5
set Interpolate output = true
set Number of grouped files = 1
end
subsection Topography
set Output to file = true
set Time between text output = 2e5
end
end
World Builder Model prm:
#### World Builder parameters which can be used by ASPECT ######## World Builder parameters which can be used by ASPECT ####
set World builder file = sub2col_v4.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-4
set Number of cheap Stokes solver steps = 100
end
end
set Output directory = sub2col_base_v68
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
set Compositional field methods = particles, particles, particles, particles, particles
set Mapped particle properties = Oceanic_Plate: initial Oceanic_Plate, Mantle: initial Mantle, Plateau_Crust: initial Plateau_Crust, Plateau_Mantle: initial Plateau_Mantle, Oceanic_Crust: initial Oceanic_Crust
end
subsection Geometry model
set Model name = box
subsection Box
set X extent = 4620e3
set Y extent = 660e3
# set Z extent = 660e3
set X repetitions = 231
set Y repetitions = 33
# set X periodic = true
end
end
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
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 = normal
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 = composite
# 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 diffusion creep = 1.5e-16 , 1.5e-16, 1.5e-16, 8.88e-15, 1.5e-16, 1e-28
set Activation energies for diffusion creep = 375.e3 , 375.e3, 375.e3, 375.e3, 375.e3, 0
set Activation volumes for diffusion creep = 3.e-6 , 3.e-6, 3.e-6, 6.e-6, 3.e-6, 0
set Prefactors for dislocation creep = 1.1e-16, 1.1e-16, 1.1e-16, 8.57e-28, 1.1e-16, 1e-50
set Stress exponents for dislocation creep = 3.5, 3.5, 3.5, 4.0, 3.5, 1
set Activation energies for dislocation creep = 530.e3, 530.e3, 530.e3, 167.e3, 530.e3, 0
set Activation volumes for dislocation creep = 14.e-6, 14.e-6, 14.e-6, 14.e-6, 14.e-6, 0
# Plasticity parameters
set Angles of internal friction = 5.7392, 5.7392, 5.7392, 5.7392, 5.7392, 0
set Cohesions = 1e6, 1e6, 1e6, 1e6, 1e6, 5e6
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, particles
subsection Visualization
set List of output variables = density, viscosity, strain rate, error indicator
set Output format = vtu
set Time between graphical output = 2e5
end
subsection Topography
set Output to file = true
set Time between text output = 5e5
end
subsection Particles
set Number of particles = 1400000
set Time between data output = 1e6
set Data output format = vtu
set List of particle properties = velocity, initial composition
set Interpolation scheme = cell average
set Update ghost particles = true
set Particle generator name = random uniform
end
end
World Builder Model wb:
{
"version":"0.4",
"cross section":[[0,330e3],[4620e3,330e3]],
"coordinate system":{"model":"cartesian"},
"features":
[
{"model":"oceanic plate", "name":"Subducting Plate", "max depth":100e3,
"coordinates":[[0,0],[0,660e3],[2600e3,660e3],[2600e3,0e3]],
"temperature models":
[
{"model":"plate model", "max depth":100e3, "spreading velocity":0.01,
"ridge coordinates":[[0e3,-1],[-100e3,2000e3]]}
],
"composition models":[
{"model":"uniform", "compositions":[4], "min depth":0, "max depth":8e3},
{"model":"uniform", "compositions":[0], "min depth":8e3,"max depth":100e3}]
},
{"model":"continental plate", "name":"Plateau_2", "coordinates":[[600e3,0],[600e3,660e3],[1000e3,660e3],[1000e3,0]],
"temperature models":[{"model":"linear", "max depth":120e3, "bottom temperature":1573}],
"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":[[3000e3,0],[3000e3,660e3],[2600e3,660e3],[2600e3,0]],
"temperature models":[{"model":"linear", "max depth":120e3, "bottom temperature":1573}],
"composition models":[
{"model":"uniform", "compositions":[2], "min depth":0, "max depth":20e3},
{"model":"uniform", "compositions":[3], "min depth":20e3,"max depth":120e3}]
},
{"model":"oceanic plate", "name":"Overriding Oceanic Plate", "max depth":100e3,
"coordinates":[[3000e3,0],[3000e3,660e3],[4620e3,660e3],[4620e3,0e3]],
"temperature models":
[
{"model":"plate model", "max depth":100e3, "spreading velocity":0.01,
"ridge coordinates":[[4620e3,-1],[4620e3, 660e3]]}
],
"composition models":[
{"model":"uniform", "compositions":[4], "min depth":0, "max depth":8e3},
{"model":"uniform", "compositions":[0], "min depth":8e3,"max depth":100e3}]
},
{"model":"mantle layer", "name":"upper mantle", "min depth":100e3, "max depth":660e3,
"coordinates":[[0,0],[0,660e3],[4620e3,660e3],[4620e3,0]],
"temperature models":[{"model":"uniform", "temperature":1573}],
"composition models":[{"model":"uniform", "compositions":[1]}]
},
{"model":"subducting plate", "name":"Slab",
"coordinates":[[2600e3,660e3],[2600e3,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":1573}],
"composition models":[
{"model":"uniform", "compositions":[4], "min distance slab top":0, "max distance slab top":8e3},
{"model":"uniform", "compositions":[0], "min distance slab top":8e3,"max distance slab top":100e3}] }
]
}
Both models have similar temperature gradients with uniform mantle temperatures of 1573 K. There is an obvious gradual viscosity increase in the first model where you can observe this easily in the viscosity field but you can’t on the other. What might be causing this?
