Implementation of a particle-in-cell method for the energy solver in 3D spherical geodynamic modeling

The thermal evolution of the Earth’s interior and its dynamic effects are the focus of Earth sciences. However, the commonly adopted grid-based temperature solver is usually prone to numerical oscillations, especially in the presence of sharp thermal gradients, such as when modeling subducting slabs and rising plumes. This phenomenon prohibits the correct representation of thermal evolution and may cause incorrect implications of geodynamic processes. After examining several approaches for removing these numerical oscillations, we show that the Lagrangian method provides an ideal way to solve this problem. In this study, we propose a particle-in-cell method as a strategy for improving the solution to the energy equation and demonstrate its effectiveness in both one-dimensional and three-dimensional thermal problems, as well as in a global spherical simulation with data assimilation. We have implemented this method in the open-source finite-element code CitcomS, which features a spherical coordinate system, distributed memory parallel computing, and data assimilation algorithms.