# Inner core convection#

*This section was contributed by Juliane Dannberg, and the model setup was
inspired by discussions with John Rudge. Additional materials and comments by
Mathilde Kervazo and Marine Lasbleis.*

This is an example of convection in the inner core of the Earth. The model is based on a spherical geometry, with a single material. Three main particularities are constitutive of this inner core dynamics modeling: it consists of a sphere where the gravity decreases linearly (to mimic self-gravitation) from the boundary to zero at the center of the inner core; the boundary conditions combine normal stress and normal velocity, and take into account the rate of phase change (melting/freezing) at the inner-outer core boundary; the material has a temperature dependent density that makes the density profile unstably stratified as temperature increases towards the center of the core. Note that we do not actually compute self-gravitation, but instead define a linear gravity profile. Since the density variations are very small, this is a good approximation.

The setup is analogous to the models described in Deguen *et al.* [2013],
and all material properties are chosen in a way so that the
equations are non-dimensional.

The required heating model and changes to the material model are implemented in the shared library inner_core_convection.cc. To compile the file, do

```
cmake -DAspect_DIR=/path/to/aspect/build/ .
make
```

In the non-dimensional form of the equations derived by Deguen *et al.* [2013],
we solve for the potential temperature
\(T = \tilde{T}-T_{\text{is}}\) (\(\tilde{T}\) is the temperature field,
\(T_{\text{is}}\) the isentropic – also called adiabatic –
temperature). This allows to solve the temperature field with simple boundary
conditions (\(T=0\)), even if the temperature of the inner core boundary evolves
with time, defined as the intersection between the isentrope and the liquidus
of the material in the outer core. The equations for inner core convection in
the approximation of no growth (equation 59 for the potential temperature) are

where \(Ra\) is the Rayleigh number and \(H\) is the ’source term’, constructed when removing the adiabatic temperature from the temperature field to obtain the potential temperature \(T\). \(H\) describes the evolution of the adiabatic temperature over time, due to secular cooling of the outer core. In spherical geometry, \(H=6\).

**Mechanical boundary.** The mechanical boundary conditions for the inner core
are tangential stress-free and continuity of the normal stress at the
inner-outer core boundary. For the non-dimensional equations, that means that
we define a “phase change number” \(\mathcal{P}\) (see Deguen *et al.* [2013]) so that the normal stress at the
boundary is \(-\mathcal{P} u_r\) with the radial velocity \(u_r\). This number
characterizes the resistance to phase change at the boundary, with
\(\mathcal{P}\rightarrow\infty\) corresponding to infinitely slow
melting/freezing (or a free slip boundary), and \(\mathcal{P}\rightarrow0\)
corresponding to instantaneous melting/freezing (or a zero normal stress,
corresponding to an open boundary).

In the weak form, this results in boundary conditions of the form of a surface integral:

with the normal vector \(\mathbf n\).

This phase change term is added to the matrix in the
inner_core_assembly.cc plugin by using a
signal (as described in Section Extending ASPECT through signals). The
signal connects the function `set_assemblers_phase_boundary`

, which is only
called once at the beginning of the model run. It creates the new assembler
`PhaseBoundaryAssembler`

for the boundary faces of the Stokes system and adds
it to the list of assemblers executed in every time step. The assembler
contains the function `phase_change_boundary_conditions`

that loops over all
faces at the model boundary, queries the value of \(\mathcal{P}\) from the
material model, and adds the surface integral given above to the matrix:

```
/*
Copyright (C) 2011 - 2023 by the authors of the ASPECT code.
This file is part of ASPECT.
ASPECT is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
ASPECT is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ASPECT; see the file LICENSE. If not see
<http://www.gnu.org/licenses/>.
*/
#include <aspect/simulator_access.h>
#include <aspect/global.h>
#include <aspect/simulator.h>
#include <aspect/simulator/assemblers/interface.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/fe/fe_values.h>
#include "inner_core_convection.cc"
namespace aspect
{
/**
* A new assembler class that implements boundary conditions for the
* normal stress and the normal velocity that take into account the
* rate of phase change (melting/freezing) at the inner-outer core
* boundary. The model is based on Deguen, Alboussiere, and Cardin
* (2013), Thermal convection in Earth's inner core with phase change
* at its boundary. GJI, 194, 1310-133.
*
* The mechanical boundary conditions for the inner core are
* tangential stress-free and continuity of the normal stress at the
* inner-outer core boundary. For the non-dimensional equations, that
* means that we can define a 'phase change number' $\mathcal{P}$ so
* that the normal stress at the boundary is $-\mathcal{P} u_r$ with
* the radial velocity $u_r$. This number characterizes the resistance
* to phase change at the boundary, with $\mathcal{P}\rightarrow\infty$
* corresponding to infinitely slow melting/freezing (free slip
* boundary), and $\mathcal{P}\rightarrow0$ corresponding to
* instantaneous melting/freezing (zero normal stress, open boundary).
*
* In the weak form, this results in boundary conditions of the form
* of a surface integral:
* $$\int_S \mathcal{P} (\mathbf u \cdot \mathbf n) (\mathbf v \cdot \mathbf n) dS,$$,
* with the normal vector $\mathbf n$.
*
* The function value of $\mathcal{P}$ is taken from the inner core
* material model.
*/
template <int dim>
class PhaseBoundaryAssembler :
public aspect::Assemblers::Interface<dim>,
public SimulatorAccess<dim>
{
public:
void
execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override
{
internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);
const Introspection<dim> &introspection = this->introspection();
const FiniteElement<dim> &fe = this->get_fe();
const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
const unsigned int n_q_points = scratch.face_finite_element_values.n_quadrature_points;
//assemble force terms for the matrix for all boundary faces
if (scratch.cell->face(scratch.face_number)->at_boundary())
{
scratch.face_finite_element_values.reinit (scratch.cell, scratch.face_number);
for (unsigned int q=0; q<n_q_points; ++q)
{
const double P = Plugins::get_plugin_as_type<const MaterialModel::InnerCore<dim>>
(this->get_material_model()).resistance_to_phase_change
.value(scratch.material_model_inputs.position[q]);
for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/)
{
if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
{
scratch.phi_u[i_stokes] = scratch.face_finite_element_values[introspection
.extractors.velocities].value(i, q);
++i_stokes;
}
++i;
}
const Tensor<1,dim> normal_vector = scratch.face_finite_element_values.normal_vector(q);
const double JxW = scratch.face_finite_element_values.JxW(q);
// boundary term: P*u*n*v*n*JxW(q)
for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
data.local_matrix(i,j) += P *
scratch.phi_u[i] *
normal_vector *
scratch.phi_u[j] *
normal_vector *
JxW;
}
}
}
};
template <int dim>
void set_assemblers_phase_boundary(const SimulatorAccess<dim> &simulator_access,
Assemblers::Manager<dim> &assemblers)
{
AssertThrow (Plugins::plugin_type_matches<const MaterialModel::InnerCore<dim>>
(simulator_access.get_material_model()),
ExcMessage ("The phase boundary assembler can only be used with the "
"material model 'inner core material'!"));
assemblers.stokes_system_on_boundary_face.push_back (std::make_unique<PhaseBoundaryAssembler<dim>>());
}
}
template <int dim>
void signal_connector (aspect::SimulatorSignals<dim> &signals)
{
signals.set_assemblers.connect (&aspect::set_assemblers_phase_boundary<dim>);
}
ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>,
signal_connector<3>)
```

Instructions for how to compile and run models with a shared library are given in Section Running benchmarks that require code.

**Governing parameters.** Analyzing
Equations (36), two
parameters determine the dynamics of convection in the inner core: the
Rayleigh number \(Ra\) and the phase change number \(\mathcal{P}\). Three main
areas can be distinguished: the stable area, the plume convection area and the
translation mode of convection area (Fig. 104). For low Rayleigh
numbers (below the critical value \(Ra_c\)), there is no convection and thermal
diffusion dominates the heat transport. However, if the inner core is
convectively unstable (\(Ra\)>\(Ra_c\)), the convection regime depends mostly on
\(\mathcal{P}\). For low \(\mathcal{P}\) (\(<29\)), the convective translation mode
dominates, where material freezes at one side of the inner core and melts at
the other side, so that the velocity field is uniform, pointing from the
freezing to the melting side. Otherwise, at high \(\mathcal{P}\) (\(>29\)),
convection takes the usual form of thermal convection with shear free boundary
and no phase change, that is the one-cell axisymmetric mode at the onset, and
chaotic plume convection for larger Rayleigh number. In this case, melting and
solidification at the ICB have only a small dynamic effect. At intermediate
values of \(\mathcal{P}\), the first unstable mode is a linear combination of the high-\(\mathcal{P}\)
convection mode and of the small-\(\mathcal{P}\) translation mode.

Changing the values of \(Ra\) and \(\mathcal{P}\) in the input file allows switching between the different regimes. The Rayleigh number can be changed by adjusting the magnitude of the gravity:

```
# The gravity has its maximum value at the boundary of inner and
# outer core, and decreases approximately linearly to zero towards
# the center of the core.
# The Rayleigh number used in the model is given by the magnitude
# of the gravity at the inner core/outer core boundary.
subsection Gravity model
set Model name = radial linear
subsection Radial linear
set Magnitude at bottom = 0.0
set Magnitude at surface = 2 # <-- Ra
end
end
```

The phase change number is implemented as part of the material model, and as a function that can depend on the spatial coordinates and/or on time:

```
subsection Material model
set Model name = inner core material
# The 'inner core material' model also contains a function that
# represents the resistance to melting/freezing at the inner core
# boundary.
# For P-->infinity, the boundary is a free slip boundary, and for
# P-->0, the boundary is an open boundary (with zero normal stress).
subsection Inner core
subsection Phase change resistance function
set Variable names = x,y,z
set Function expression = 1e-2 # <-- P
end
end
end
```

Fig. 105 shows examples of the three regimes with \(Ra=3000, \mathcal{P}=1000\) (plume convection), \(Ra=10^5, \mathcal{P}=0.01\) (translation), \(Ra=10, \mathcal{P}=30\) (no convection).

**Mesh refinement.** The temperature is set to 0 at the outer boundary and a
large temperature gradient can develop at the boundary layer, especially for
the translation regime. The adaptive mesh refinement allows it to resolve this
layer at the inner core boundary. Another solution is to apply a specific
initial refinement, based on the boundary layer thickness scaling law
\(\delta \propto Ra^{-0.236}\), and to refine specifically the uppermost part of
the inner core.

In order to have a mesh that is much finer at the outer boundary than in the center of the domain, this expression for the mesh refinement subsection can be used in the input file:

```
subsection Mesh refinement
set Initial global refinement = 4 #this may be more expensive, and should be run on a cluster.
set Initial adaptive refinement = 1
set Strategy = minimum refinement function
set Time steps between mesh refinement = 0
subsection Minimum refinement function
set Variable names = depth, phi, theta
set Function expression = if(depth>0.1,if(depth>0.2,2,5),6)
end
end
```

**Scaling laws.** In addition, Deguen *et al.* [2013]
give scaling laws for the velocities in each regime derived from linear
stability analysis of perfect translation, and show how numerical results
compare to them. In the regimes of low \(\mathcal{P}\), translation will start
at a critical ratio of Rayleigh number and phase change number
\(\frac{Ra}{\mathcal{P}}=\frac{175}{2}\) with steady-state translation
velocities being zero below this threshold and tending to
\(v_0=\frac{175}{2}\sqrt{\frac{6}{5}\frac{Ra}{\mathcal{P}}}\) going towards
higher values of \(\frac{Ra}{\mathcal{P}}\). In the same way, translation
velocities will decrease from \(v_0\) with increasing \(\mathcal{P}\), with
translation transitioning to plume convection at \(\mathcal{P}\sim29\). Both
trends are shown in Fig. 106 and can be compared to Figure 8
and 9 in Deguen *et al.* [2013].

Deguen, R., T. Alboussière, and P. Cardin. 2013. “Thermal
Convection in Earth’s Inner Core with Phase Change at Its
Boundary.” *Geophysical Journal International* 194 (3): 1310–34.
https://doi.org/10.1093/gji/ggt202.