Mandyoc: A finite element code to simulate thermochemical convection in parallel

Mandyoc is a 2-D finite element code written in C dedicated to simulating thermochemical convection in the interior of terrestrial planets. Different linear and non-linear rheologies can be adopted, appropriately simulating the strain and stress pattern in the Earth’s crust and mantle, both in extensional and collisional tectonics. Additionally, the code allows variations of boundary condition for the velocity field in space and time, simulating different pulses of tectonism in the same numerical scenario.


Summary
Mandyoc is a 2-D finite element code written in C dedicated to simulating thermochemical convection in the interior of terrestrial planets. Different linear and non-linear rheologies can be adopted, appropriately simulating the strain and stress pattern in the Earth's crust and mantle, both in extensional and collisional tectonics. Additionally, the code allows variations of boundary condition for the velocity field in space and time, simulating different pulses of tectonism in the same numerical scenario.

Statement of need
Mandyoc, the acronym for MANtle DYnamics simulatOr Code, is designed to simulate Stokes flow type thermochemical convection taking different compositional layers into account, and it is also appropriate to simulate Earth's lithospheric dynamics on geological timescales. Similar codes are available for geodynamic problems, like ASPECT (Bangerth et al., 2021b(Bangerth et al., , 2021aHeister et al., 2017;Kronbichler et al., 2012) and Underworld (Beucher et al., 2019;Moresi et al., 2007). Therefore, Mandyoc is an alternative to preexistent softwares.
One advantage of Mandyoc is the possibility to create scenarios with velocity boundary conditions variable in space and time, allowing the user to simulate different tectonic pulses in the same model run. Additionally, the current version incorporates surface processes, imposing rates of erosion and sedimentation on the top of the free surface (Silva & Sacek, 2022). Recently, other thermomechanical codes available for the scientific community incorporated the interaction with surface processes (e.g., Beucher & Huismans, 2020;Neuharth et al., 2022) taking into account the simulation of fluvial and hillslope processes. In the present version of Mandyoc, only predefined erosion/sedimentation rate (variable in space and time) is possible.
Previous versions of the code were used to study the evolution of continental margins, showing the interaction of the continental lithosphere with the asthenospheric mantle (Sacek, 2017;Salazar-Mora & Sacek, 2021).

Mathematics
Mandyoc solves the equations for conservation of mass, momentum and energy using the Finite Element Method assuming the extended Boussinesq approximation, respectively: u i is the component i of the velocity field, T is temperature, t is time, κ is the thermal diffusivity, H is the volumetric heat production, c p is the specific heat capacity, g is gravity, ρ is the effective rock density dependent on temperature and composition, ρ 0 is the reference rock density at temperature T 0 , α is the coefficient of thermal expansion, P is the dynamic pressure, η is the effective viscosity, and δ ij is the Kronecker delta.
The code is fully parallelized using the Portable, Extensible Toolkit for Scientific Computation (PETSc) (Balay et al., 1997(Balay et al., , 2021a(Balay et al., , 2021b. The present version of the code can simulate thermochemical convection using different rheological formulations: Newtonian flow, non-linear viscous flow or visco-plastic deformation. For example, the lithosphere can be simulated as a combination of different visco-plastic layers in which the effective viscosity depends on a nonlinear power law viscous rheology and a plastic yield criterion, like the Drucker-Prager criterion. Additionally, strain softening is implemented to facilitate the localization of strain in the plastic regime during, for example, lithospheric stretching. The composition and strain history is tracked by particles present in the interior of the finite element. The exchange of particles among the subdomains of the model is efficiently parallelized in PETSc using DMSwarm (May & Knepley, 2017).
The free surface of the Earth can be simulated and is numerically stabilized using the Free Surface Stabilization Algorithm (Kaus et al., 2010). Surface processes of erosion and sedimentation can also be incorporated in the thermo-mechanical model, allowing the coupling between the Earth's interior dynamics and the processes occurring at its surface. Complex boundary conditions for the velocity field, variable both in space and time, can be adopted by the user to simulate different episodes of tectonism. Different benchmarks are available in the repository and can be reproduced by the user, including thermochemical convection (Van Keken et al., 1997) and plume-lithosphere interaction (Crameri et al., 2012).
As an example application of Mandyoc, Figure 1 presents snapshots of one numerical scenario of lithospheric stretching imposing a divergent flow direction, resulting in rifting and break-up. In this example, the upper crust, lower crust, lithospheric mantle and asthenosphere present different rheology and density, resulting in faulting mainly in the upper crust and part of the lithospheric mantle. Additionally, deformations in the lower crust and at the base of the lithospheric mantle are accommodated by ductile creep flow. This example can be reproduced from the repository.