# Cookbook for a global-scale 2D box mantle convection model
# with melt migration.
# In this file we will go through all of the steps that are
# required for adding melting and melt transport to a mantle
# convection simulation.

# For models with melt migration, there is a nonlinear coupling between
# the Stokes system, the temperature, and the advection equation for the
# porosity (several material properties, such as the viscosities and the
# permeability depend nonlinearly on the porosity; and changes in temperature
# determine how much material is melting or freezing).
# Because of that, we use a nonlinear solver scheme ('iterated Advection and Stokes')
# that iterates between all of these equations, and we have to set its
# solver tolerance and the maximum number of iterations separately from
# the linear solver parameters.
set Nonlinear solver scheme                = iterated Advection and Stokes
set Max nonlinear iterations               = 20
set Nonlinear solver tolerance             = 1e-5

# In addition, melting and freezing normally happens on a much faster
# time scale than the flow of melt, so we want to decouple the advection
# of melt (and temperature) and the melting process itself. To do that,
# we use the operator splitting scheme, and define that for every
# advection time step, we want to do at least 10 reaction time steps.
# If these time steps would be larger than 10,000 years, we will do
# more reaction time steps (so that reaction time step size never exceeds
# 10,000 years). Here, we also specify the Stokes linear solver tolerance
# and maximum number of cheap Stokes solver steps to improve the nonlinear
# convergence behavior.
set Use operator splitting                     = true

subsection Solver parameters
  subsection Operator splitting parameters
    set Reaction time step                     = 1e4
    set Reaction time steps per advection step = 10
  end

  subsection Stokes solver parameters
    set Linear solver tolerance = 1e-8
    set Number of cheap Stokes solver steps = 100
  end
end

subsection Melt settings
  # In addition, we now also specify in the model settings that we want to
  # run a model with melt transport.
  set Include melt transport                  = true
end

##################### Settings for melt transport ########################

# In models with melt transport, we always need a compositional field with
# the name 'porosity'. Only the field with that name will be advected with
# the melt velocity, all other compositional fields will continue to work
# as before. Material models will typically query for the field with the
# name porosity to compute all melt material properties.
# In addition, the 'melt global' material model also requires a field with the
# name 'peridotite'. This field is used to track how much material has been
# molten at each point of the model, so it tracks the information how the
# composition of the rock changes due to partial melting events (sometimes
# also called depletion). This is important, because usually less melt is
# generated for a given temperature and pressure if the rock has undergone
# melting before. Typically, material properties like the density are also
# different for more or less depleted material.
subsection Compositional fields
  set Number of fields = 2
  set Names of fields = porosity, peridotite
end

##################### Initial conditions ########################

# Now that our model uses compositional fields, we also need initial
# conditions for the composition.
# We use the function plugin to set both fields to zero at the beginning
# of the model run.
subsection Initial composition model
  set Model name = function

  subsection Function
    set Function expression = 0; 0
    set Variable names      = x,y
  end
end

##################### Boundary conditions ########################

# We again choose the initial composition as boundary condition
# for all compositional fields.
subsection Boundary composition model
  set List of model names = initial composition
end

# Models with melt transport also need an additional boundary condition:
# the gradient of the fluid pressure at the model boundaries. This boundary
# condition indirectly also prescribes boundary conditions for the melt velocity,
# as the melt velocity is related to the fluid pressure gradient via Darcy's law.
# If we choose the fluid pressure gradient = solid density * gravity, melt will
# flow in and out of the model (even if the solid can not flow out) according to
# the dynamic fluid pressure in the model. Conversely, if we choose the
# fluid pressure gradient = fluid density * gravity, melt will flow in or out
# with the same velocity as the solid (so for a closed boundary, no melt will
# flow in or out). This is what we choose as our boundary condition here.
subsection Boundary fluid pressure model
  set Plugin name = density

  subsection Density
    set Density formulation = fluid density
  end
end

##################### Material properties ########################

# In addition to the material properties for the solid rock,
# we also have to specify properties for the melt.
subsection Material model
  set Model name = melt global

  subsection Melt global
    # First we describe the parameters for the solid, in the same way
    # we did in the model without melt transport
    set Thermal conductivity              = 4.7
    set Reference solid density           = 3400
    set Thermal expansion coefficient     = 2e-5
    set Reference shear viscosity         = 5e21
    set Thermal viscosity exponent        = 7
    set Reference temperature             = 1600
    set Solid compressibility             = 4.2e-12

    # The melt usually has a different (lower) density than the solid.
    set Reference melt density            = 3000

    # The permeability describes how well the pores of a porous material
    # are connected (and hence how fast melt can flow through the rock).
    # It is computed as the product of the reference value given here
    # and the porosity cubed. This means that the lower the porosity is
    # the more difficult it is for the melt to flow.
    set Reference permeability            = 1e-8

    # The bulk viscosity describes the resistance of the rock to dilation
    # and compaction. Melt can only flow into a region that had no melt
    # before if the matrix of the solid rock expands, so this parameter
    # also limits how fast melt can flow upwards.
    # The bulk viscosity is computed as the reference value given here times
    # a term that scales with one over the porosity. This means that for zero
    # porosity, the rock can not dilate/compact any more, which is the same
    # behavior that we have for solid mantle convection.
    set Reference bulk viscosity          = 1e19

    # In dependence of how much melt is present, we also weaken the shear
    # viscosity: The more melt is present, the weaker the rock gets.
    # This scaling is exponential, following the relation
    # viscosity ~ exp(-alpha * phi),
    # where alpha is the parameter given here, and phi is the porosity.
    set Exponential melt weakening factor = 10

    # In the same way the shear viscosity is reduced with increasing temperature,
    # we also prescribe the temperature-dependence of the bulk viscosity.
    set Thermal bulk viscosity exponent   = 7

    # Analogous to the compressibility of the solid rock, we also define a
    # comressibility for the melt (which is generally higher than for the solid).
    # As we do not want our compressibility to depend on depth, we set the
    # pressure derivative to zero.
    set Melt compressibility              = 1.25e-11
    set Melt bulk modulus derivative      = 0.0

    # Finally, we prescribe the viscosity of the melt, which is used in Darcy's
    # law. The lower this viscosity, the faster melt can flow.
    set Reference melt viscosity          = 1

    # change the density contrast of depleted material (in kg/m^3)
    set Depletion density change          = -200.0

    # How much melt has been generated and subsequently extracted from a particular
    # volume of rock (how 'depleted' that volume of rock is) usually changes the
    # solidus. The more the material has been molten already, the less melt will be
    # generated afterwards for the same pressure and temperature conditions. We
    # model this using a simplified, linear relationship, saying that to melt 100%
    # of the rock the temperature has to be 200 K higher than to melt it initially.
    set Depletion solidus change          = 200

    # We also have to determine how fast melting and freezing should happen.
    # Here, we choose a time scale of 10,000 years, which is a relatively long time
    # (or in other words, slow melting rate), but because this is a global model
    # and the time steps are big, it should be sufficient.
    set Melting time scale for operator splitting = 1e4
  end
end

##################### Mesh refinement #########################

# For the model with melt migration, we use adaptive refinement.
# We make use of two different refinement criteria: we set a minimum of 5 global
# refinements everywhere in the model (which is the same resolution as for the
# model without melt), and we refine in regions where melt is present, to be
# precise, everywhere where the porosity is bigger than 1e-6.
# We adapt the mesh every 5 time steps.
subsection Mesh refinement
  set Coarsening fraction                      = 0.05
  set Refinement fraction                      = 0.8
  set Initial adaptive refinement              = 2
  set Initial global refinement                = 5
  set Strategy                                 = composition threshold, minimum refinement function, boundary
  set Time steps between mesh refinement       = 5

  # minimum of 5 global refinements
  subsection Minimum refinement function
    set Function expression = 5
  end

  # refine where the porosity is bigger than 1e-6
  subsection Composition threshold
    set Compositional field thresholds = 1e-6,1.0
  end

  # refine at top and bottom boundary
  subsection Boundary
    set Boundary refinement indicators = top, bottom
  end
end

##################### Postprocessing ########################

# In addition to the visualization output, we select a number
# of postprocessors that allow us to compute some statistics
# about the output (to see how much the model without and the
# model with melt migration differ), and in particular we use
# the "depth average" postprocessor that will allow us to plot
# depth-averaged model quantities over time.
subsection Postprocess
  set List of postprocessors = visualization, composition statistics, velocity statistics, temperature statistics, depth average

  # For the model with melt migration, also add a visualization
  # postprocessor that computes the material properties relevant
  # to migration (permeability, viscosity of the melt, etc.).

  subsection Visualization
    set List of output variables      = material properties, nonadiabatic temperature, strain rate, melt material properties

    subsection Material properties
      set List of material properties = density, viscosity, thermal expansivity
    end

    subsection Melt material properties
      set List of properties = fluid density, permeability, fluid viscosity, compaction viscosity
    end
  end
end
