Paul James Tackley

A desirable characteristic of mantle convection models is the ability to determine surface topography evolution over time on a global scale. This capability is increasingly important due to growing interest in coupling planetary geodynamics with climate, landscape, habitats, and biological evolution systems. A common way of achieving this in numerical geodynamic models with a fixed Eulerian grid is through the implementation of a free-surface boundary condition using the sticky air method in conjunction with an appropriate method of tracking the free surface. Although existing methods for tracking the interface between the air and rock layers are available, they often struggle to provide high-resolution results on a global scale without incurring significant computational costs. We propose a method for representing surfaces directly using Lagrangian markers that can track the location of a free surface in 2D and 3D models and test it using the finite-volume mantle convection code StagYY. This approach offers a direct, high-resolution representation of the surface without the need for a large number of high-cost volumetric markers throughout the model domain. Benchmarks demonstrate the effectiveness of this method compared to the commonly-used marker-in-cell method. The direct representation of the surface enables additional features such as the direct tracking of sea levels over time, potential for coupling with surface process models on the global scale, and enables the implementation of alternative discretisations of the Stokes equations to improve Stokes solver accuracy near the free surface boundary.

Accurately and efficiently modelling topographic evolution is a key challenge in geodynamic modelling, which requires the solution of the Stokes equations with free surface boundary conditions. While finite difference methods on staggered grids, as used in geodynamic modelling codes such as StagYY, I3ELVIS and LaMEM, offer strong computational performance and compatibility with multigrid solvers, the use of fixed Eulerian grids complicates the implementation of realistic, deformable free surfaces. Two existing methods are available to model free surface boundary conditions in such codes: the commonly used sticky-air method, which suffers from limitations relating to high viscosity contrasts, and the “staircase” method, which improves upon the sticky air method by imposing free surface boundary conditions at cell boundaries. To address the limitations of existing methods of implementing free surface boundary conditions, this study investigates an alternative variational discretisation of the Stokes equations that uses volume fractions to represent a smooth surface within a fixed Eulerian grid, allowing the imposition of accurate free surface boundary conditions while allowing it to bypass the limitations of existing free surface discretisation methods. The variational Stokes method is demonstrated to be an accurate and computationally efficient alternative to existing methods. It reproduces results comparable to existing methods while reducing computational cost and enabling broader applications, including non-zero surface tractions, complex surface loading, and compatibility with 3D spherical geometries.

A desirable characteristic of mantle convection models is the ability to determine surface topography evolution over time on a global scale. This capability is increasingly important due to growing interest in coupling planetary geodynamics with climate, landscape, habitats, and biological evolution systems. A common way of achieving this in numerical geodynamic models with a fixed Eulerian grid is through the implementation of a free-surface boundary condition using the sticky air method in conjunction with an appropriate method of tracking the free surface. Although existing methods for tracking the interface between the air and rock layers are available, they often struggle to provide high-resolution results on a global scale without incurring significant computational costs. We propose a method for representing surfaces directly using Lagrangian markers that can track the location of a free surface in 2D and 3D models and test it using the finite-volume mantle convection code StagYY. This approach offers a direct, high-resolution representation of the surface without the need for a large number of high-cost tracers throughout the model domain. Benchmarks demonstrate the effectiveness of this method compared to the commonly-used marker-in-cell method. The direct representation of the surface enables additional features such as the direct tracking of sea levels over time, potential for coupling with surface process models on the global scale, and enables the implementation of alternative discretisations of the Stokes equations to improve Stokes solver accuracy near the free surface boundary.

A desirable characteristic of mantle convection models is the ability to determine surface topography evolution over time on a global scale. This capability is increasingly important due to growing interest in coupling planetary geodynamics with climate, landscape, habitats, and biological evolution systems. A common way of achieving this in numerical geodynamic models with a fixed Eulerian grid is through the implementation of a free-surface boundary condition using the sticky air method in conjunction with an appropriate method of tracking the free surface. Although existing methods for tracking the interface between the air and rock layers are available, they often struggle to provide high-resolution results on a global scale without incurring significant computational costs. We propose a method for representing surfaces directly using Lagrangian markers that can track the location of a free surface in 2D and 3D models and test it using the finite-volume mantle convection code StagYY. This approach offers a direct, high-resolution representation of the surface without the need for a large number of high-cost tracers throughout the model domain. Benchmarks demonstrate the effectiveness of this method compared to the commonly-used marker-in-cell method. The direct representation of the surface enables additional features such as the direct tracking of sea levels over time, potential for coupling with surface process models on the global scale, and enables the implementation of alternative discretisations of the Stokes equations to improve Stokes solver accuracy near the free surface boundary.

Accurately and efficiently modelling topographic evolution is a key challenge in geodynamic modelling, which requires the solution of the Stokes equations with free surface boundary conditions. While finite difference methods on staggered grids, as used in geodynamic modelling codes such as StagYY, I3ELVIS and LaMEM, offer strong computational performance and compatibility with multigrid solvers, the use of fixed Eulerian grids complicates the implementation of realistic, deformable free surfaces. Two existing methods are available to model free surface boundary conditions in StagYY: the commonly used sticky-air method, which suffers from limitations relating to high viscosity contrasts, and the "staircase" method, which improves upon the sticky air method by imposing free surface boundary conditions at cell boundaries. To address the limitations of existing methods of implementing free surface boundary conditions, this study investigates an alternative variational discretisation of the Stokes equations that uses volume fractions to represent a smooth surface within a fixed Eulerian grid, allowing the imposition of accurate free surface boundary conditions while allowing it to bypass the limitations of existing free surface discretisation methods. The variational Stokes method is demonstrated to be an accurate and computationally efficient alternative to existing methods. It reproduces results comparable to existing methods while reducing computational cost and enabling broader applications, including non-zero surface tractions, complex surface loading, and compatibility with 3D spherical geometries.

Accurately and efficiently modelling topographic evolution is a key challenge in geodynamic modelling, which requires the solution of the Stokes equations with free surface boundary conditions. While finite difference methods on staggered grids, as used in geodynamic modelling codes such as StagYY, I3ELVIS and LaMEM, offer strong computational performance and compatibility with multigrid solvers, the use of fixed Eulerian grids complicates the implementation of realistic, deformable free surfaces. Two existing methods are available to model free surface boundary conditions in StagYY: the commonly used sticky-air method, which suffers from limitations relating to high viscosity contrasts, and the "staircase" method, which improves upon the sticky air method by imposing free surface boundary conditions at cell boundaries. To address the limitations of existing methods of implementing free surface boundary conditions, this study investigates an alternative variational discretisation of the Stokes equations that uses volume fractions to represent a smooth surface within a fixed Eulerian grid, allowing the imposition of accurate free surface boundary conditions while allowing it to bypass the limitations of existing free surface discretisation methods. The variational Stokes method is demonstrated to be an accurate and computationally efficient alternative to existing methods. It reproduces results comparable to existing methods while reducing computational cost and enabling broader applications, including non-zero surface tractions, complex surface loading, and compatibility with 3D spherical geometries.

Radiogenic heating plays a crucial role in shaping a planet’s evolution and dynamics. On Earth, ~50% of surface heat loss originates from the decay of three long-lived, heat-producing elements (HPEs): potassium, thorium, and uranium. These elements are strongly lithophile and preferentially concentrate in the silicate mantle of planets. However, a recent study by Luo et al. (Science Advances, 2024) suggests that under the high-pressure, high-temperature conditions of core formation in large rocky planets (so-called super-Earths), these HPEs may become siderophile, partitioning preferentially into the iron core. The presence of HPEs in the mantles of super-Earths plays a crucial role in their internal dynamics. A feedback loop between internal heating, temperature, and viscosity regulates mantle temperature, adjusting viscosity to the value needed to facilitate convective loss of the radiogenic heat (Tackley et al., Icarus 2013). However, if these sources of radiogenic heat partition into the core, mantle convection in super-Earths becomes dominated by heat flowing from the core rather than by a mix of internal heating and cooling from above (as in Earth). Using 1D, parameterized mantle evolution models, Luo et al. (Science Advances, 2024) show that this shift leads to a sharp rise in core-mantle boundary (CMB) temperatures and an increase in total CMB heat flow, with significant implications for volcanism and magnetic field generation.In this study, we perform mantle convection simulations using the StagYY code (Tackley, PEPI 2008), extending the models of Tackley et al. (Icarus, 2013) to include HPEs in the core, as suggested by Luo et al. (Science Advances, 2024). Our models are run in a 2D spherical annulus geometry and allow for melting at all mantle depths. We test different planetary masses, from 1 to 10 Earth masses, as well as different post-perovskite rheologies, (upper- and lower-bound, following Tackley et al. 2013, and interstitial rheology following Karato 2011), two tectonic regimes (stagnant and mobile-lid), and three mantle-to-core partitioning ratios of HPEs (0.1, 1, and 10). This work contributes to the growing understanding of the interior dynamics of super-Earths, and their implications on surface and atmospheric conditions, the presence of a magnetic field, and habitability potential.

Giant impacts were common in the early evolution of the Solar System. Such an impact has been suggested to have affected the rotation of Venus and possibly its thermal evolution, causing long-lived volcanic activity.Here, we explore a range of possible giant impacts using smoothed particle hydrodynamics (SPH). We analyse the post-impact rotation and debris disc masses to identify scenarios that can reproduce Venus’ present-day characteristics. We model post-impact interior dynamics evolution using the StagYY convection code, by transferring the thermal field obtained through SPH simulations into a 2D spherical annulus geometry. We account for the low viscosity of molten volumes of the mantle above 35% melt fraction by using an effective "eddy" thermal conductivity of 1010 W/(m.K) following Lourenço et al. (2020), and the heat flux is parametrized following Abe (1993, 1997). The evolution of the core uses a 1D parameterized model to track the temperature profile and the growth of the inner core.We observe that a wide range of impact scenarios are consistent with Venus’ current rotation for both head-on collisions on a non-rotating Venus and oblique, hit-and-run impacts on a rotating Venus. Collisions that match consistent rotation rates typically produce minimal debris discs residing within Venus’ synchronous orbit (Bussmann et al., 2025). We select these favourable scenarios to model their long-term interior evolution.In the simulations, giant impacts expectedly produce surface magma oceans. Their relative depths vary between different simulations depending on impact properties: from a shallow melt layer in the order of 100 km thick to a fully molten mantle from surface to core-mantle boundary for the most energetic impacts (high impactor mass and velocity). If the surface is able to radiate heat to space efficiently, the magma ocean cools down quickly and first reaches the rheological transition (35% melt fraction) in a few 100-1000 yrs. Full solidification (0% melt fraction) can take longer because of the effects of the impact on the deep interior.Indeed, as highlighted by Marchi et al. (2023), giant impacts also deposit a considerable amount of energy in the upper layers of the core, which translates into temperatures reaching up to 104 K. This causes the base of the mantle to fully melt. The resulting liquid layers (in the core and the mantle) convect vigorously and cool the core rapidly (~104 years). The very hot mantle melt is buoyant and rises toward the surface through the solid mantle on timescales of 104-105 yrs. Plumes formed in such a way persist until the excess of heat is extracted from the core and the lower mantle reaches the rheological transition. Solidification of the surface can be delayed by plumes, but models indicate that a fully solid state is reached in a few 1-10 Myr.After a few hundred million years, the thermal evolution of a Venus-like planet that experienced a giant impact becomes similar to that of cases devoid of impacts. The characteristics of the convection regime in both cases do not substantially differ at present-day (after 4.5 Gyr).

Venus and Earth exemplify two divergent evolutionary pathways of rocky planets despite their similar sizes, bulk densities and orbital distances. Today, Earth maintains a thin N2- and O2-dominated atmosphere, regulated by efficient volatile and crustal recycling, while Venus hosts a CO2-dominated, 92 bar atmosphere that results in extreme surface temperatures. The origin of Venus’ massive atmosphere is likely a combination of processes including early magma ocean degassing, impact delivery and long-term mantle outgassing, with possibly one or more catastrophic outgassing events (Gillmann et al., 2022). The efficiency and composition of long-term mantle outgassing is still poorly constrained for Venus. Previous studies have shown that pressure is the primary control on volatile solubility, regulating whether degassing occurs at all and the elemental sequence of volatility e.g. C-bearing species degass predominantely at higher pressures, while the redox state of the mantle controls primary the chemical speciation of outgassed volatiles e.g. wheather CH4 or CO or CO2 is degassed (Gaillard and Scaillet, 2014, Gaillard et al. 2020). In contrast to Ortenzi et al. (2020), who have focused on the redox-dependence of volatile partitioning and gas speciation for various planetary masses, this study investigates the pressure-dependence of volcanic outgassing for varying redox states and its implications for Venus’ atmospheric evolution. First, Atmodeller (Bower et al., 2025) is used to calculate the pressure-dependent volatile solubility and gas speciation for erupted lava. Then, the findings are linked to 2D mantle convection models of Venus with varying surface temperatures and rheological parameters. These models use the numerical mantle convection code StagYY (Tackley, 2008) including partial melting with prescribed intrusion efficiencies to quantify melt production, lithospheric thickness as well as the tectonic regime. By explicitly accounting for volatile degassing as a function of surface pressure, our study aims to constrain total outgassing and the role of pressure-buffered volatile retention in shaping Venus’ atmosphere.