# 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