Long-term simulation of Martian Polar Layered Deposits: introducing the Planetary Evolution Model

Introduction:
Understanding the long-term evolution of the Martian climate and its reservoirs is a key challenge in planetary science. Mars experiences large obliquity variations (0°–60°), which are known to drive its climate [1]. In particular, the North Polar Layered Deposits (NPLD) is believed to preserve a valuable stratigraphic archive of past environmental changes over the last few to tens of millions of years [2,3]. Yet, a major question remains: linking the orbital forcings with the layers evolution under climate dynamics [4]. To investigate these processes at geologic timescales, we developed the Planetary Evolution Model (PEM), a new modeling framework to bridge the gap between short-term General Circulation Models (GCM) and long-term climate dynamics.
In this contribution, we present the conceptual and physical basis of the PEM and illustrate its first applications focused on the formation and evolution of the Martian NPLD.

The Planetary Evolution Model (PEM):
Traditional Mars GCM convincingly simulate atmospheric processes but are limited to short periods (typically decades) due to computational costs. 1D climate or stratigraphic models handle long durations but often neglect key dynamical couplings, limiting their applicability. Moreover, all these models typically rely on paramatrized fluxes and prescribed reservoirs (glaciers, subsurface ice) whose initialization can lead to unrealistic outcomes. Hence, GCM are unsuitable for simulating orbital-scale climate variations. The PEM addresses this issue by focusing on long-term changes while bypassing sub-year variability. It is based on asynchronous coupling with a GCM, here we use the Mars Planetary Climate Model (PCM) [3].
The core idea of the PEM is to compute tendencies from two consecutive PCM years and then to extrapolate them to evolve climate-relevant reservoirs until stopping criteria are met (e.g. surface pressure change, ice area loss/gain or orbital shift). At that point, the PEM halts and, with the updated climate state, it re-runs the PCM to obtain new tendencies. This cycle repeats until the desired simulation length is reached. The PEM can be used in two ways: (i) to determine steady-state configurations (e.g. realistic distribution of glaciers) for given orbital parameters, which can then be used as initial states in other simulations; and (ii) to perform transient simulations over full orbital cycles.
Key physical processes in the PEM include:
    • CO₂/H₂O surface ice evolution: local tendencies are computed from interannual minima of the ice in the PCM (perennial ice). The PEM assumes fast atmospheric equilibration for water with only transfers between sublimating and condensing reservoirs.
    • Glacier flow: a statistical sub-grid slope parameterization accounts for North–South orientation effects [6], reproducing slope-dependent layering.
    • Subsurface H₂O ice: ice table depth is dynamically adjusted based on thermal diffusion, pressure and surface humidity, following Norbert Schorghofer's work [7].
    • Soil properties: the PEM uses a multilayer subsurface (regolith, breccia, bedrock) whose thermal properties adapt to pressure and pore-filling ice. The yearly-averaged surface and subsurface temperatures are updated according to surface ice accumulation/depletion.
    • CO₂/H₂O adsorption/desorption: adsorbed species are maintained in equilibrium with the atmosphere.
Furthermore, the PEM performs several subsequent adaptations to make the aforementioned processes internally consistent and to simulate the co-evolution of reservoirs and atmosphere. For instance, the PEM adjusts the yearly-averaged surface pressure based on CO₂ mass balance. This necessitates to scale the volume mixing ratios accordingly and to correct the CO₂ tendencies obtained from the PCM.

Layered Deposits:
The NPLD exhibit a complex stratigraphic structure made of layers with varying composition (H₂O ice, dust and possibly tempoeary CO₂ ice) and geometry (thickness, unconformities). The alternation of dusty and icy layers has long been hypothesized to record obliquity-driven climate cycles. High obliquities (>30°) typically lead to ice loss, while low obliquities (<20°) enable accumulation provided that a source is available [8,9].
Additionally, CO₂ and H₂O sublimation/condensation thresholds interact to create diverse physical siturations. For instance, water may act as a lag layer above CO₂ glaciers, as observed and modeled by Buhler et al. [10] for the South Polar Layered Deposits (SPLD).
To reproduce such features, the PEM includes a novel dynamic layer-tracking model. Each deposition or ablation event modifies the local stratigraphy by adding new layers with appropriate dust/ice content or generating dust lag deposits upon ice sublimation. 
As a first scientific application, we used the PEM to simulate the growth and evolution of the NPLD under an orbital forcing scenario. The goal is to test whether the PEM can reproduce observed stratigraphic features and evaluate the role of orbital cycles in shaping the NPLD.

Perspectives:
The PEM is a new-generation tool to study the evolution of Mars’ climate and volatile reservoirs over orbital timescales and understand its present-day geomorphological imprints. Our first applications to the NPLD demonstrate the PEM potential to resolve long-term polar stratigraphy and test hypotheses on past Martian environments, bringing us closer to decoding the Mars’ deep climate memory.
Our ongoing work aims to improve the PEM models to explore atmospheric collapse/inflation scenarios and to fully couple the PEM with an hydrological model (lakes, river systems) to study early Mars.

References:
[1] Laskar J. et al. (2004) Icarus, 170(2), 343–364.
[2] Phillips R. J. et al. (2011) Science, 332, 838.
[3] Levrard B. et al. (2007) J. Geophys. Res., 112, E06012.
[4] Smith I. B. et al. (2020) Planetary and Space Science, 184, 104841.
[5] Forget F. et al. (1999) J. Geophys. Res., 104(E10), 24155–24175.
[6] Lange L. et al. (2023) Journal of Geophysical Research: Planets, 128, e2023JE007915.
[7] Schorghofer N. (2007) Nature, 449, 192-194.
[8] Hvidberg C. S. et al. (2012) Icarus, 221(1), 405-419.
[9] Vos E. et al. (2022) J. Geophys. Res., 127(3), e2021JE007115.
[10] Buhler P. B. et al. (2020) Nature Astronomy, 4, 364-371.