Mantle convection

Mantle convection¶

This code is taken from a chapter in The FEniCS Book, which is available for free online. That chapter is in turn an implementation of a model setup from van Keken et al (1997). I took the thermomechanical parts and removed the chemistry. This paper and the code from the FEniCS book all use a non-dimensional form of the equations, which I've adopted here.

In this notebook, we'll see the Stokes equations again, but we'll couple them to the evolution of temperature through an advection-diffusion equation. The extensive quantity is not the temperature itself but the internal energy density $E = \rho c_p T$ where $\rho$, $c_p$ are the mass density and the specific heat at constant pressure. The flux of heat is $$F = \rho c_p T u - k\nabla T.$$ We'll assume there are no heat sources but, but the real physics would include sources from decay of radioactive elements, strain heating, and chemical reactions. The variational form of the heat equation is: $$\int_\Omega\left\{\partial_t(\rho c_p T)\,\phi - (\rho c_p Tu - k\nabla T)\cdot\nabla\phi - Q\phi\right\}dx = 0$$ for all test functions $\phi$. I've written this in such a way that it still makes sense when the density and heat capacity aren't constant. We'll assume they are (mostly) constant in the following. The variational form of the momentum balance equation is $$\int_\Omega\left\{2\mu\varepsilon(u): \dot\varepsilon(v) - p\nabla\cdot v - q\nabla\cdot u - \rho g\cdot v\right\}dx = 0$$ for all velocity and pressure test functions $v$, $q$.

The characteristic that makes this problem interesting is that we assume the flow is just barely compressible. This is called the Boussinesq approximation). In practice, it means we do two apparently contradictory things. First, we assume that the flow is incompressible: $$\nabla\cdot u = 0.$$ Second, we assume that the density that we use on the right-hand side of the momentum balance equation actually is compressible, and that the degree of compression is linear in the temperature: $$\rho = \rho_0(1 - \beta(T - T_0))$$ where $\rho_0$ is a reference density and $\beta$ is thermal expansion coefficient. These assumptions aren't consistent with each other. We're also not consistent in using this everywhere. For example, we only use this expression on the right-hand side of the momentum balance equation, but we'll still take the internal energy density to be equal to $\rho_0c_pT$ in the energy equation. In other words, the material is compressible as far as momentum balance is concerned but not energy balance. I find that uncomfortable. But I would find writing a full compressible solver to be blinding agony. Anyway it lets us make pretty pictures.

I won't go into a full non-dimensionalization of the problem, which is detailed in the van Keken paper. The most important dimensionless numbers for this system are the Prandtl number and the Rayleigh number. The Prandtl number quantifies the ratio of momentum diffusivity to thermal diffusivity: $$\text{Pr} = \frac{\mu / \rho}{k / \rho c_p} = \frac{\nu}{\kappa}$$ where $\nu$ and $\kappa$ are respectively the kinematic viscosity and thermal diffusivity. For mantle flows, the Prandtl number is on the order of $10^{25}$. We can get a sense for a characteristic velocity scale from the momentum balance equation. First, we can take a factor of $\rho_0g$ into the pressure. If we take $U$, $L$, and $\Delta T$ to be the characteristic velocity, length, and temperature difference scales, non-dimensionalizing the momentum balance equation gives $$\frac{\mu U}{L^2} \approx \rho_0g\beta\Delta T$$ which then implies that $U = \rho_0g\beta L^2\Delta T / \mu$. Now we can compute the thermal Péclet number $$\text{Pe} = \frac{UL}{\kappa}.$$ When the Péclet number is large, we're in the convection-dominated regime; when it's small, we're in the diffusion-dominated regime. Taking our estimate of velocity scale from before, we get $$\text{Pe} = \frac{\rho_0 g \beta L^3 \Delta T}{\kappa\mu}.$$ This is identical to the Rayleigh number, which we'll write as $\text{Ra}$. Lord Rayleigh famously used linear stability analysis to show that a fluid that is heated from below and cooled from above is unstable when the Rayleigh number exceeds 600 or so. Here we'll take the Rayleigh number to be equal to $10^6$ as in the van Keken paper. It's a worthwhile exercise in numerical analysis to see what it takes to make a solver go for higher values of the Rayleigh number.

Here we're using a rectangular domain. Again, the spatial scales have been non-dimensionalized.

The initial condition for temperature involves loads of annoyingly long expressions, but it's all in Appendix A of the van Keken paper.

The important point is to make the fluid hotter below and on one side and cooler above and on the other side, shown in the plot below.

Next we need to make some function spaces for the fluid velocity and pressure. Note how the degree of the velocity space is one higher than that of the pressure -- we're using the Taylor-Hood element again.

Once we've created a function in the mixed space, we can then pull out the two parts with the split method.

The code below creates the variational form of the Stokes equations $$\nabla\cdot\tau - \nabla p + \text{Ra}\;T\; g = 0$$ with temperature-dependent gravitational forcing.

We can use the .sub method to pull parts out of mixed spaces, which we need in order to create the right boundary conditions.

A bit of magic in order to tell the linear solver that the Stokes equations have a null space we need to project out.

Here I'm creating solver object since we'll need to repeatedly solve the same system many times.

Here we're setting up the temperature solver, which includes both convection and diffusion. I've set the timestep based on the mesh size and the maximum speed that we just found above from the velocity solution.

And the timestepping loop. Note that the final time is on a non-dimensional scale again, in physical time it works out to be on the order of a hundred million years. Here we're using an operator splitting approach. We first update the temperature, and then we compute a new velocity and pressure. The splitting error goes like $\mathscr{O}(\delta t)$ as the timestep is reduced. So if we use a first-order integration scheme like backward Euler then the splitting error is asymptotically the same as that of the discretization itself. We could get $\mathscr{O}(\delta t^3)$ or higher convergence by using, say, the Radau-IIA method, but the total error would be dominated by splitting, so there would be little point in trying harder.

The movie below shows the temperature evolution. The rising and sinking plumes of hot and cold fluid correspond to the mantle plumes that produce surface volcanism.

What next¶

There are a few features in the van Keken paper and otherwise that I didn't attempt here. The original van Keken paper adds chemistry. Solving the associated species transport equation with as little numerical diffusion as possible is hard. You could then make the fluid buoyancy depend on chemical composition and add temperature-dependent chemical reactions.

Real mantle fluid also has a temperature-dependent viscosity. When I tried to add this using the splitting scheme above, the solver quickly diverges. Getting that to work might require a fully coupled time-integration scheme, which I did try. Rather than split up the temperature and the velocity/pressure solves, a fully coupled scheme would solve for all three fields at once. If you could make this work, it would open up the possibility of using higher-order methods like Radau-IIA. Whether it would go or not seemed to depend on the machine and the day of the week. I'll explore that more in a future post.