swh:1:snp:75cdaf5164207cb3d00f07a3da10a0250b29d03b
Tip revision: 7e867c1f3ae18b4869f32f42694e83e30e7e5051 authored by Timo Heister on 11 January 2021, 19:20:39 UTC
Merge pull request #3946 from gassmoeller/fix_changelog
Merge pull request #3946 from gassmoeller/fix_changelog
Tip revision: 7e867c1
mid_ocean_ridge.prm
# 2d model for melting and melt transport at a mid-ocean ridge.
# As the flow at mid-ocean ridges can be assumed to be roughly symmetric
# with respect to the ridge axis in the center, we only model one half
# of the ridge.
set Dimension = 2
set Adiabatic surface temperature = 1570
# Because our model includes melt transport, it is nonlinear and we have to
# use an iterative solver scheme, that iterates between solving the temperature
# composition and Stokes equations.
set Nonlinear solver scheme = iterated Advection and Stokes
set Output directory = output_mid_ocean_ridge
# The end time of the simulation. We want to run the model until it reaches
# steady state, which is after approximately 6 million years.
set End time = 8e6
##################### Melting and freezing ########################
# Because the model includes reactions that might be on a faster time scale
# than the time step of the model (melting and the freezing of melt), we use
# the operator splitting scheme.
set Use operator splitting = true
subsection Solver parameters
subsection Operator splitting parameters
# We choose the size of the reaction time step as 200 years, small enough
# so that it can accurately model melting and freezing.
set Reaction time step = 2e2
# Additionally, we always want to do at least 10 operator splitting time
# steps in each model time step, to accurately compute the reactions.
set Reaction time steps per advection step = 10
end
# Because this model includes strong localized viscosity contrasts we
# increase the robustness of the solver at the cost of memory consumption.
subsection Stokes solver parameters
set GMRES solver restart length = 200
end
end
# We use the melt simple material model that includes melting and freezing of
# melt for an average mantle composition that is characteristic for a mid-ocean
# ridge setting, and mainly use its default parameters.
# In particular, we have to define how fast melting and freezing should be.
# We assume that both reactions happen on a time scale of 200 years (or a rate
# of 5e-3/year), which should be substantially shorter than the time step size,
# so that the melt fraction will always be close to equilibrium.
# As the model includes melting and freezing, we do not have to extract any melt.
# We also slightly reduce the dependence of shear viscosity on the melt
# fraction, because the model is relatively coarse and would otherwise develop
# velocity oscillations in the region of highest melt content.
subsection Material model
set Model name = melt simple
subsection Melt simple
set Reference permeability = 1e-7
set Melt extraction depth = 0.0
set Freezing rate = 0.005
set Melting time scale for operator splitting = 2e2
set Exponential melt weakening factor = 20
end
end
##################### Model geometry ########################
# Our model geometry is a box of 105x70 km. This guarantees that inflowing
# material is solid, and will start to melt within the model domain.
subsection Geometry model
set Model name = box
subsection Box
set X extent = 105000
set Y extent = 70000
# To keep the aspect ratio of our elements close to one, we chose the
# coarse mesh is such a way that it has more elements in X than in Y
# direction, in the same ratio as the aspect ratio of the model.
set X repetitions = 3
set Y repetitions = 2
end
end
# The gravity is constant and points downward.
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 10.0
end
end
##################### Velocity ########################
# To model the divergent velocitiy field of a mid-ocean ridge, we prescribe
# the plate velocity (pointing away from the ridge) at the top boundary.
# We use a closed boundary with free slip conditions as the left boundary, which
# marks the ridge axis and also acts as a center line for our model, so that
# material can not cross this boundary.
# We prescribe the velocity at the top boundary using a function:
# At the ridge axis, the velocity is zero, at a distance of 10 km from the ridge
# axis or more, the rigid plate uniformly moves away from the ridge with a constant
# speed, and close to the ridge we interpolate between these two conditions.
subsection Boundary velocity model
set Prescribed velocity boundary indicators = top:function
set Tangential velocity boundary indicators = left
subsection Function
# We choose a half-spreading rate of u0=3cm/yr.
set Function constants = u0=0.03, x0=10000
set Variable names = x,z
set Function expression = if(x<x0,(1-(x/x0-1)*(x/x0-1))*u0,u0); 0
end
end
# We prescribe the lithostatic pressure as a boundary traction on
# the bottom and right side of the model, so that material can flow in and out
# according to the flow induced by the moving plate.
subsection Boundary traction model
set Prescribed traction boundary indicators = right:initial lithostatic pressure, bottom:initial lithostatic pressure
subsection Initial lithostatic pressure
# We calculate the pressure profile at the right model boundary.
set Representative point = 105000, 70000
end
end
##################### Temperature ########################
# As initial temperature, we choose an adiabatic profile with boundary layers at the
# top and the bottom. We make the top boundary layer very old (100 Ma) so that in the
# beginning, the material is still solid and the porosity is zero.
subsection Initial temperature model
set Model name = adiabatic
subsection Adiabatic
set Age top boundary layer = 1e8
set Age bottom boundary layer = 1e5
set Amplitude = 0
subsection Function
set Function expression = 0;0
end
end
end
# We choose a constant temperature at the top and bottom of the domain.
# In particular, the bottom boundary temperature will control the steady state
# temperature profile, and we choose a potential temperature of 1570 K (1300 °C).
subsection Boundary temperature model
set Fixed temperature boundary indicators = top, bottom
set List of model names = box
subsection Box
set Top temperature = 293
set Bottom temperature = 1570
end
end
# We want to include latent heat of melting and freezing in the model, as it is
# an essential process for the temperature evolution, but as the model domain
# is quite small and this is a simple model, we do not include any other sources
# of heat.
subsection Heating model
set List of model names = latent heat
end
# Melt moves with a different velocity than the solid, and transports energy,
# so we include this process in the model.
subsection Melt settings
# We want to solve the McKenzie Equations to track the flow of melt.
set Include melt transport = true
set Heat advection by melt = true
end
##################### Composition ########################
# We need two compositional fields: The porosity field to track the motion of
# melt, and the peridotite field to track the depletion of material, which is
# changed by the melting and freezing reactions.
subsection Compositional fields
set Number of fields = 2
set Names of fields = porosity, peridotite
end
# We set both fields to zero at the start of the model: The material is solid
# (zero porosity) and has the average mantle composition (zero depletion).
subsection Initial composition model
set Model name = function
subsection Function
set Function expression = 0; 0
set Variable names = x,y
end
end
# The boundary conditions (which are relevant for material inflowing at the
# bottom boundary of the model) are the same as the initial conditions.
subsection Boundary composition model
set Fixed composition boundary indicators = bottom
set List of model names = initial composition
end
##################### Mesh refinement #########################
# We use adaptive mesh refinement to increase the resolution in regions where
# melt is present, and otherwise use a uniform grid.
subsection Mesh refinement
set Coarsening fraction = 0.5
set Refinement fraction = 0.5
# A refinement level of 5 (4 global + 1 adaptive refinements) corresponds to
# a cell size of approximately 1 km.
set Initial adaptive refinement = 1
set Initial global refinement = 4
set Strategy = minimum refinement function, composition threshold
set Time steps between mesh refinement = 5
subsection Minimum refinement function
set Coordinate system = cartesian
set Function expression = 4
set Variable names = x,y
end
# We use a very small refinement threshold for the porosity to make sure that
# all cells where the two-phase flow equations are solved (melt cells) have
# the higher resolution.
subsection Composition threshold
set Compositional field thresholds = 1e-6, 1.0
end
end
##################### Postprocessing ########################
subsection Postprocess
set List of postprocessors = visualization, composition statistics, velocity statistics
# We mainly want to look at material properties of the solid and the melt.
subsection Visualization
set List of output variables = material properties, melt material properties, melt fraction
subsection Material properties
set List of material properties = density, viscosity
end
# To see in which cells melt transport is modelled, it can be useful to look
# at the property 'is melt cell', so we include it in the output. In addition,
# we always visualize the compaction pressure 'p_c' if this postprocessor is
# used.
subsection Melt material properties
set List of properties = compaction viscosity, permeability, fluid density, is melt cell
end
set Time between graphical output = 0
end
end
# We write a checkpoint every 100 time steps, so that we are able to restart
# the computation from that point.
subsection Checkpointing
set Steps between checkpoint = 100
end