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

(36)#\[\begin{split} \begin{aligned} \nabla \cdot \sigma &= -Ra T \mathbf g, \\ \nabla \cdot \mathbf u &= 0, \\ \left(\frac{\partial T}{\partial t} + \mathbf u\cdot\nabla T\right) - \nabla^2 T &= H, \end{aligned}\end{split}\]

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:

\[\int_S \mathcal{P} (\mathbf u \cdot \mathbf n) (\mathbf v \cdot \mathbf n) \text{d}S,\]

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 the signals mechanism). 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 - 2022 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:

      virtual
      void
      execute (internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
               internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
      {
        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. 99). 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.

Screenshot

Fig. 99 Stability diagram for convection in a sphere with phase change at its outer boundary. The stability curves for the first unstable mode (\(l=1\)) and the translation are obtained from Deguen et al. [2013]. Each dot (no convection) and triangle (blue: translation, yellow: plume convection) is one model run done with ASPECT. The higher \(Ra\) and \(\mathcal{P}\) are, the more refinement is required (see text).#

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 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-->inifinity, 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. 100 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).

Screenshot

Fig. 100 Convection regimes in the inner core for different values of \(Ra\) and \(\mathcal{P}\). From left to right: no convection, translation, plume convection; the 2D slices in the top row use the default temperature scale for all panels, while in the second row an adaptive scale is used. The bottom row features slightly different model parameters (that are still in the same regime as the models in the respective panels above) and also shows the velocity as arrows.#

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. 101 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.