Numerical solution of a non-linear conservation law applicable to the interior dynamics of partially molten planets

TL;DR

This paper develops a robust numerical scheme to solve a non-linear diffusion equation modeling heat transport in partially molten planets, capturing large entropy gradient variations and applicable to planetary interior dynamics.

Contribution

It introduces a stable finite volume method for non-linear diffusion equations with large dynamic ranges, specifically tailored for planetary magma ocean modeling.

Findings

Successfully captures large entropy gradient variations (~12 orders of magnitude).

Provides a flexible, energy-conserving numerical framework applicable to planetary science.

Enables extended precision calculations for complex non-linear diffusion problems.

Abstract

The energy balance of a partially molten rocky planet can be expressed as a non-linear diffusion equation using mixing length theory to quantify heat transport by both convection and mixing of the melt and solid phases. In this formulation the effective or eddy diffusivity depends on the entropy gradient, $\partial S/\partial r$, as well as entropy. First we present a simplified model with semi-analytical solutions, highlighting the large dynamic range of $\partial S/\partial r$, around 12 orders of magnitude, for physically-relevant parameters. It also elucidates the thermal structure of a magma ocean during the earliest stage of crystal formation. This motivates the development of a simple, stable numerical scheme able to capture the large dynamic range of $\partial S/\partial r$ and provide a flexible and robust method for time-integrating the energy equation. We then consider a full model including energy fluxes associated with convection, mixing, gravitational separation, and conduction that all depend on the thermophysical properties of the melt and solid phases. This model is discretised and evolved by applying the finite volume method (FVM), allowing for extended precision calculations and using $\partial S/\partial r$ as the solution variable. The FVM is well-suited to this problem since it is naturally energy conserving, flexible, and intuitive to incorporate arbitrary non-linear fluxes that rely on lookup data. Special attention is given to the numerically challenging scenario in which crystals first form in the centre of a magma ocean. Our computational framework is immediately applicable to modelling high melt fraction phenomena in Earth and planetary science research. Furthermore, it provides a template for solving similar non-linear diffusion equations arising in other disciplines, particularly for non-linear functional forms of the diffusion coefficient.