Hi all,
I’m trying to study the influence of plume-lithosphere interaction on rift formation. But when I add the mantle plume as a thermal anomaly( ΔT=150°C) in the model, the viscosity distribution of the initial model will be uneven in the horizontal direction.
Model without mantle plume
Model with mantle plume
I guess there is something wrong with the setting of mantle plume, but I don’t know what the specific problem is. Here is my prm file.
#### Global parameters
set Dimension = 2
set Start time = 0
set End time = 20e6
set Use years in output instead of seconds = true
set Nonlinear solver scheme = single Advection, iterated Stokes
set Nonlinear solver tolerance = 1e-4
set Max nonlinear iterations = 10 # Normally set to 50 or 100
set CFL number = 0.5
set Maximum time step = 50e3 # Normally set to 20 Kyr
set Use conduction timestep = true
set Output directory = test_convection_wz50
set Timing output frequency = 1
set Pressure normalization = no
#subsection Termination criteria
# set Termination criteria = end step
# set End step = 100
#end
#subsection Formulation
# set Formulation = Boussinesq approximation
#end
# Model geometry (20 x 20km)
subsection Geometry model
set Model name = box
subsection Box
set X extent = 800e3
set Y extent = 400e3
set X repetitions = 8
set Y repetitions = 4
end
end
# Mesh (2.5 x 2.5km)
subsection Mesh refinement
set Initial global refinement = 3
set Initial adaptive refinement = 0
set Time steps between mesh refinement = 0
# set Strategy = viscosity
# set Strategy = minimum refinement function
# subsection Minimum refinement function
# set Coordinate system = cartesian
# set Variable names = x,y
# set Function expression = if ( x>=700.e3 && x<=900.e3, 4, 3)
# end
set Run postprocessors on initial refinement = false
end
# Quadratic elements for compostion, temperature, and Stokes velocity
subsection Discretization
set Composition polynomial degree = 2
set Stokes velocity polynomial degree = 2
set Temperature polynomial degree = 2
end
# Solver parameters
subsection Solver parameters
subsection Stokes solver parameters
set Stokes solver type = block AMG
set Linear solver tolerance = 1e-5
set Number of cheap Stokes solver steps = 0
end
# subsection Newton solver parameters
# set Max Newton line search iterations = 5
# set Max pre-Newton nonlinear iterations = 10
# set Maximum linear Stokes solver tolerance = 1e-1
# set Nonlinear Newton solver switch tolerance = 1e-4
# set SPD safety factor = 0.9
# set Stabilization preconditioner = SPD
# set Stabilization velocity block = SPD
# set Use Newton failsafe = false
# set Use Newton residual scaling method = false
# set Use Eisenstat Walker method for Picard iterations = true
# end
end
# A normal free surface advection may become unstable when significant topography develops
subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface
subsection Free surface
set Surface velocity projection = normal
end
end
####### boundary condition #######
subsection Boundary velocity model
set Tangential velocity boundary indicators = left, right, bottom
end
# Number and names of compositional fields
subsection Compositional fields
set Number of fields = 6
set Names of fields = Upper_crust, Lower_crust, Lithosphere_mantle, Mantle, Weak_zone, Mantle_plume
end
# Spatial domain of different compositional fields
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = \
if(y>=380.e3, 1, 0); \
if(y<380.e3 && y>=360.e3, 1, 0); \
if((y<360.e3 && y>=280.e3 && x>=0 && x<375.e3) || (y<360.e3 && y>=280.e3 && x>425.e3 && x<=800.e3) , 1, 0); \
if(y<280.e3 && y>=0 && sqrt((x-400.e3)^2 + y^2)>50.e3, 1, 0); \
if(x>=375.e3 && x<=425.e3 && y>=280.e3 && y<360.e3, 1, 0); \
if( sqrt((x-400.e3)^2 + y^2)<=50.e3 && y>=0, 1, 0)
end
end
# Composition: fixed on bottom, free on sides and top
subsection Boundary composition model
set Fixed composition boundary indicators = bottom
set List of model names = initial composition
end
# Temperature boundary conditions
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set Model name = box
subsection Box
set Bottom temperature = 1720
set Left temperature = 0
set Right temperature = 0
set Top temperature = 273
end
end
# Initial temperature field
# Surface Temperature - upper crust (ts1) = 273 K
# Surface Heat Flow - upper crust (qs1) = 0.055 W/m^2
# Heat Production - upper crust (A1) = 1.00e-6 W/m^3;
# Heat Production - lower crust (A2) = 0.25e-6 W/m^3;
# Heat Production - mantle_Lith (A3) = 0.00e-6 W/m^3;
# Thermal Conductivity - crust = 2.5 ; lith_m =3.0 (W/(m K));
# To satisfy these constraints, the following values are required:
# Surface Temperature - Lower crust (ts2) = 633 K
# - Mantle_Lithosphere (ts3) = 893 K
# - mantle (ts4) = 1593K
# Surface Heat Flow - Lower crust (qs2) = 0.035 W/m^2;
# - Mantle_Lithosphere (qs3) = 0.02625 W/m^2;
subsection Initial temperature model
set List of model names = function
subsection Function
set Variable names = x,y
set Function constants = h=400.e3,ts1=273,ts2=633,ts3=893,Tp=1300,Cp=1250,a=3.0e-5,g=9.81, \
A1=1.00e-6,A2=0.25e-6,A3=0.0, \
k1=2.5,k2=2.5,k3=3.0, \
qs1=0.055,qs2=0.035,qs3=0.02625
set Function expression = if( (h-y)<=20.e3, \
ts1 + (qs1/k1)*(h-y) - (A1*(h-y)*(h-y))/(2.0*k1), \
if( (h-y)>20.e3 && (h-y)<=40.e3, \
ts2 + (qs2/k2)*(h-y-20.e3) - (A2*(h-y-20.e3)*(h-y-20.e3))/(2.0*k2), \
if( (h-y)>40.e3 && (h-y)<=120.e3, \
ts3 + (qs3/k3)*(h-y-40.e3) - (A3*(h-y-40.e3)*(h-y-40.e3))/(2.0*k3) ,\
273+Tp*exp(a*g*(h-y)/Cp)+ if( sqrt((x-400.e3)^2 + y^2) <= 50.e3 && y>=0, 150, 0) ) ) ) ;
end
end
# Constant internal heat production values (W/m^3) for background material
# and compositional fields.
subsection Heating model
set List of model names = compositional heating
subsection Compositional heating
set Compositional heating values = 0.0, 1.0e-6, 0.25e-6, 0.0, 0.0, 0.0, 0.0
end
end
# Material model
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
# Reference temperature
set Reference temperature = 273
set Minimum strain rate = 1.e-20
set Reference strain rate = 1.e-16
# Limit the viscosity with minimum and maximum values
set Minimum viscosity = 1e18
set Maximum viscosity = 1e26
set Define thermal conductivities = true
# Upper_c Lower_c Lith_m Mantle Weak_z Mantle_p
set Thermal conductivities = 3.3, 2.5, 2.5, 3.0, 3.3, 3.0, 3.3
set Heat capacities = 1250, 750, 750, 1250, 1250, 1250, 1250
# Δp=150
set Densities = 3300, 2800, 2900, 3300, 3300, 3300, 3150
set Thermal expansivities = 3.0e-5, 2.5e-5, 2.5e-5, 3.0e-5, 3.0e-5, 3.0e-5, 3.0e-5
# Harmonic viscosity averaging
set Viscosity averaging scheme = harmonic
set Viscous flow law = composite
# Dislocation creep parameters for Upper_c Lower_c Lith_m Mantle Weak_z Mantle_p
set Prefactors for dislocation creep = 1.10e-16, 1.00e-15, 2.06e-23, 1.10e-16, 1.10e-16, 1.10e-16, 1.10e-16
set Stress exponents for dislocation creep = 3.5, 2.0, 3.2, 3.5, 3.5, 3.5, 3.5
set Activation energies for dislocation creep = 530.e3, 167.e3, 238.e3, 530.e3, 530.e3, 530.e3, 530.e3
set Activation volumes for dislocation creep = 20.e-6, 0., 0., 17.e-6, 20.e-6, 20.e-6, 20.e-6
set Prefactors for diffusion creep = 2.46e-16, 1.00e-15, 2.06e-23, 2.46e-16, 2.46e-16, 2.46e-16, 2.46e-16
set Stress exponents for diffusion creep = 1.0, 2.0, 3.2, 1.0, 1.0, 1.0, 1.0
set Activation energies for diffusion creep = 375.e3, 167.e3, 238.e3, 375.e3, 375.e3, 375.e3, 375.e3
set Activation volumes for diffusion creep = 10.e-6, 0., 0., 10.e-6, 10.e-6, 10.e-6, 10.e-6
set Grain size = 1e-3
set Grain size exponents for diffusion creep = 3.0, 1.0, 1.0, 1.0, 3.0, 1.0, 3.0
set Angles of internal friction = 20., 20., 20., 20., 20., 3., 20.
set Cohesions = 20.e6, 20.e6, 20.e6, 20.e6, 20.e6, 2.e6, 20.e6
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 = basic statistics, composition statistics, heat flux densities, heat flux statistics, mass flux statistics, matrix statistics, pressure statistics, temperature statistics, topography, velocity statistics, visualization
subsection Visualization
set List of output variables = density, heat flux map, named additional outputs, strain rate, viscosity, principal stress, stress
set Time between graphical output = 100.e3
set Interpolate output = false
set Output format = vtu
end
subsection Topography
set Output to file = true
end
end