# Sublithospheric convection beneath an oceanic plate with a grain size dependent rheology#

*This section was contributed by Juliane Dannberg.*

The input file for this model can be found at cookbooks/grain_size_ridge/grain_size_ridge.prm

Note that this model may take a significant time to run due to the high resolution and large number of time steps. We recommend to either run it on a workstation or cluster or to reduce the model resolution.

This model features the spreading of an oceanic plate away from a mid-ocean ridge, with the ridge being in the top left corner of the model. Its vertical extent is 410 km and and its horizontal extent is 4000 km. Since the plate moves with a velocity of 5 cm/yr, this leads to a maximum plate age of 80 Myr. Material flows in from the bottom, and then leaves the domain at the right side. The initial temperature is the adiabatic temperature of the mantle with an added cold boundary layer becoming thicker with increasing distance from the ridge (increasing plate age). The mineral grain size starts from a constant value of 5 mm and then evolves over time. Grains grow faster at higher temperature, and the grain size is reduced by deformation in the dislocation creep regime.

## Grain size evolution.#

To model the evolution of the mineral grain size and to take this grain size
into account in our rheological model, we use the grain size material model,
where these equations are implemented. This model implements a composite
diffusion–dislocation rheology, with diffusion creep depending on the
evolving grain size, and deformation by dislocation creep reducing the grain
size. For more details on this material model and how to use it for global
mantle convection models, see Dannberg *et al.* [2017].
We have to specify the parameters for both creep mechanisms and for the grain
size evolution.

```
subsection Material model
set Model name = grain size
subsection Grain size model
# Diffusion creep has a prefactor, a temperature dependence defined by the
# activation energy, a pressure dependence defined by the activation volume
# and a grain size dependence. We take the activation energy from Hirth &
# Kohlstedt, 2003, and pick the prefactor and activation volume so that we
# get a reasonable upper mantle viscosity profile and a balance between
# diffusion and dislocation creep.
set Diffusion creep prefactor = 5e-15
set Diffusion activation energy = 3.75e5
set Diffusion activation volume = 4e-6
set Diffusion creep grain size exponent = 3
# Dislocation creep has a prefactor, a temperature dependence defined by the
# activation energy, a pressure dependence defined by the activation volume
# and a strain-rate dependence defined by the exponent. We take the activation
# energy and volume from Hirth & Kohlstedt, 2003, and pick the prefactor so
# that we get a reasonable upper mantle viscosity and a balance between
# diffusion and dislocation creep.
set Dislocation creep prefactor = 1e-15
set Dislocation activation energy = 5.3e5
set Dislocation activation volume = 1.4e-5
set Dislocation creep exponent = 3.5
set Dislocation viscosity iteration number = 10000
# Grain size is reduced when work is being done by dislocation creep.
# Here, 10% of this work goes into reducing the grain size rather than
# into shear heating.
# Grain size increases with a rate controlled by the grain growth
# constant and a temperature-depndence defined by the activation
# energy. By setting the activation volume to zero, we disable the
# pressure-dependence of grain growth.
# The parameter values are taken from Dannberg et al., 2017, G-cubed,
# with the parameters being based on Faul and Jackson, 2007.
set Work fraction for boundary area change = 0.1
set Grain growth rate constant = 1.92e-10
set Grain growth activation energy = 4e5
set Grain growth activation volume = 0
# Viscosity is cut off at a minimum value of 10^16 Pa s
# and a maximum value of 10^23 Pa s.
set Maximum viscosity = 1e23
set Minimum viscosity = 1e16
set Maximum temperature dependence of viscosity = 1e8
end
end
```

Many of these parameters, in particular the prefactors and activation volumes, have large uncertainties. The values used here should therefore only be seen as example values. It often makes sense to choose their values based on the viscosity profile one wants to use. These viscosity profiles are usually constrained by observations, such as the geoid or mineral physics data, and therefore provide additional constraints to the values we get from experiments.

Note that diffusion creep depends on the third power of the grain size as
specified by the `Diffusion creep grain size exponent`

parameter, and that
the stress-dependence of dislocation creep is controlled by the
`set Dislocation creep exponent`

parameter.

## Tracking the grain size on particles.#

To track this evolving mineral grain size on particles, we have to have a
compositional field with the name `grain_size`

, we have to specify that
we want to advect this field using the particle method, and we have to select
`grain_size`

as the particle property to be mapped to this compositional field:

```
# Our only compositional field is the grain size, and we use particles
# to track it. In order for the grain size to not just be tracked, but
# to evolve over time depending on the temperature and deformation in
# the model, we select ``grain_size'' as the mapped particle property
# for the grain_size compositional field.
subsection Compositional fields
set Number of fields = 1
set Names of fields = grain_size
set Compositional field methods = particles
set Mapped particle properties = grain_size:grain_size
end
```

In addition, we have to specify a number of different parameters controlling how to initialize and advect these particles, how to interpolate the computed values back to the mesh and how to control the load balancing. We here choose to add and remove particles in such a way that there are always at least 40 and not more than 640 particles in each cell. This means that a cell with 160 particles can either be refined once (resulting in about 40 particles in each new cell) or coarsened once (resulting in about 640 in the new bigger cell) without having to add or remove particles. We also use the bilinear least squares interpolation scheme to interpolate from the particles to the grid, which is more accurate than schemes based on, for example, the cell average. To make sure this scheme does not generate any over- or under- shoots in the grain size, we switch on the linear least squares limiter. To use this scheme, we also have to update the ghost particles, because the interpolation scheme requires knowledge about particles in neighboring cells. Finally, we use a second Order Runge Kutta integrator to advect the particles.

```
subsection Postprocess
set List of postprocessors = visualization, composition statistics, velocity statistics, mass flux statistics, particles
# We use particles to advect the grain size.
subsection Particles
set Interpolation scheme = bilinear least squares
set Number of particles = 500000
set Time between data output = 1e6
set List of particle properties = grain size
set Load balancing strategy = remove and add particles
set Minimum particles per cell = 40
set Maximum particles per cell = 640
set Integration scheme = rk2
set Update ghost particles = true
subsection Interpolator
subsection Bilinear least squares
set Use linear least squares limiter = true
end
end
end
end
```

## Model evolution.#

We run the model for 200 million years. In the first few tens of millions of years there are no large changes in the temperature structure, except that the thermal boundary layer increases slightly in thickness. In addition, the grain size evolves from its initially constant value, with the grain size being reduced due to deformation. This leads to a band of particularly low grain size in the lowermost part of the oceanic plate. Grain size is also reduced in the asthenesphere, but to a lesser degree.

At 58.8 Myr, the first downwelling from the base of the lithosphere starts to form. Due to the strong temperature-dependence of the viscosity, only the lowermost part of the plate participates in this downwelling, which is therefore very thin compared to the plate. Because of the strong deformation within the downwelling (which reduces the grain size) and its lower temperature (which hinders grain growth), the grain size in these convective instabilities is smaller compared to the surrounding mantle. Over time, the downwelling moves towards the left and younger plate ages, and additional downwellings start to develop until the model reaches a quasi-steady state.

Note that for simplicity, our model uses a constant thermal conductivity of 4.7 W/m/K, whereas the thermal conductivity is likely lower in the lithosphere and is also temperature-dependent. Therefore, our lithosphere is thicker (and downwellings occur earlier) than on Earth.