# Glacial isostatic adjustment reveals Mars's interior viscosity structure

A. Broquet, A. -C. Plesa, V. Klemann, B. C. Root, A. Genova, M. A.

Wieczorek, M. Knapmeyer, J. C. Andrews-Hanna, D. Breuer

► To cite this version:

A. Broquet, A. -C. Plesa, V. Klemann, B. C. Root, A. Genova, et al. Glacial isostatic adjustment reveals Mars's interior viscosity structure. Nature, 2025, 639, pp.109-113. 10.1038/s41586-024-08565-9. insu-05080916

HAL Id: insu-05080916

https://insu.hal.science/insu-05080916v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Glacial isostatic adjustment reveals Mars's interior viscosity structure

https://doi.org/10.1038/s41586-024-08565-9

Received: 20 August 2024

Accepted: 18 December 2024

Published online: 26 February 2025

Open access

<figure><img src="image_2.png" /></figure>

Check for updates

A. Broquet¹², A.-C. Plesa<fnref n="1" />, V. Klemann², B. C. Root³, A. Genova⁴, M. A. Wieczorek⁵, M. Knapmeyer<fnref n="1" />, J. C. Andrews-Hanna⁶ & D. Breuer<fnref n="1" />

Investigating glacial isostatic adjustment has been the standard method to decipher Earth's interior viscosity structure<fnref n="1" />,², but such an approach has been rarely applied to other planets because of a lack of observational data³,⁴. The north polar cap of Mars is the only millions-of-years-old surface feature that can induce measurable surface deformation on this planet, thereby holding clues to its present-day internal viscosity structure⁵,⁶. Here we investigate the emplacement of this ice cap by combining thermal evolution models⁷, viscoelastic deformation calculations⁸ and radar observations⁶. We show that downward motion of the northern regions is ongoing and can be constrained by analyses of the time-variable gravity field⁹ and NASA's InSight seismic moment rate¹⁰. Only models with present-day high viscosities (2–6 × 10²² Pa s for depths greater than 500 km), strong mantle depletion in radiogenic elements (more than 90%) and thick average crusts (thicker than 40 km) are consistent with the negligible flexure beneath the polar cap seen by radars. The northern lithosphere must deform at less than 0.13 mm per year and have a seismic efficiency less than 0.3 to satisfy gravity and seismic constraints, respectively. Our models show that the north polar cap formed over the last 1.7–12.0 Myr and that glacial isostatic adjustment can be further constrained by future gravity recovery missions to Mars¹¹,¹².

The response of a planet to loading is intimately linked to its interior structure<fnref n="1" />,². On Earth, studying the time-variable response of the lithosphere to the growth and decay of large ice sheets, or glacial isostatic adjustment, has been the standard approach to constrain our planet's upper mantle viscosity structure<fnref n="1" />. However, such models have rarely been applied to other terrestrial bodies due to a lack of observational data³,⁴. Mars harbours two geologically young (less than 100 million years old) and large (more than 1,000 km across) polar ice caps¹³, which represent the only millions-of-years-old surface features that can induce measurable viscoelastic deformation on this planet. In the absence of in situ heat flow measurements, the analysis of these deformations is one of the few methods that give access to the present-day thermal state and interior structure of Mars⁵⁻⁷.

Orbital radar sounders have mapped the Martian ice caps⁵,¹⁴ and the lack of clear measurable lithospheric flexure beneath the north polar deposits (Fig. 1) was used to constrain the planet's interior to be cold with a thick and stiff elastic lithosphere⁵,⁶. Whereas the results of these studies are used as constraints on the present-day thermal state of Mars's interior¹⁵,¹⁶, geodynamic evolution models⁷,¹⁷,¹⁸ struggle to explain the thick and cold lithosphere inferred at the north pole (Extended Data Fig. 1) and, at the same time, the planet's young volcanism¹⁹,²⁰ and ongoing plume activity²¹. This indicates that the assumptions of elastic flexure in previous models should be revisited to account for the interior's transient viscoelastic response.

Most previous analyses have assumed the polar deformation to be at equilibrium, which is only valid if the time elapsed since the polar cap formed is greater than the time required for viscous adjustments⁵,⁶. Geologic observations and global climate models suggest the north polar cap is only a few million years old, but the exact age remains uncertain²²,²³. Owing to this young age, calculations suggested that viscoelastic relaxation could result in estimating thinner elastic lithospheres and higher heat flows beneath the north polar cap⁵,⁶. However, these models assumed a constant mantle viscosity, did not account for the ice loading history and only considered a single wavelength to represent the load. As a result, limited insights into the possible effects of viscoelastic relaxation were provided. Here we investigate glacial isostatic adjustment on Mars in light of newly acquired constraints on the planet's interior structure from NASA's InSight mission⁷,¹⁵,²⁴⁻²⁶. Thermal evolution models⁷ are used to parameterize the viscoelastic interior structure and are further constrained by the long-term loading history of the north polar cap and radar observations. The north polar deformations are limited based on the observed planetary time-variable gravity field⁹ and the InSight-derived seismic moment rate¹⁰, allowing to provide tight constraints on Mars's interior viscosity structure.

## Mars's interior structure from geodynamic models

Without plate tectonics on Mars, incompatible radioactive heat sources, which were sequestered in the crust by magmatic processes,

<footnote n="1">¹Institute for Planetary Research, German Aerospace Center, DLR, Berlin, Germany. ²GFZ Helmholtz Centre for Geosciences, Potsdam, Germany. ³Delft University of Technology, Delft, The Netherlands. ⁴Sapienza University of Rome, Rome, Italy. ⁵Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, Paris, France. ⁶University of Arizona, Tucson, AZ, USA.
⁷e-mail: adrien.broquet@dlr.de</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1 | Deformation beneath the north polar cap as seen by radars.**

a,b, Context image of the north pole (a) and interpolated deformed basement assuming a real dielectric constant of 3.0 based on 78 SHARAD and 213 MARSIS radar measurements (b) (Methods). c, SHARAD depth-corrected radargram (orbit 51924) showing negligible downward flexure. The radargram apparent</figcaption></figure>

depth was corrected by assigning real permittivities of 1.0 and 3.0 above and below the detected ground surface. The radar track (i,j) is shown in **a**, **b** and **c**. We note that because of its higher frequency, SHARAD does not penetrate through the sand-rich basal unit beneath the centre of the north polar cap⁵.

are unlikely to have been reincorporated into the mantle in any substantial amount²⁷. This naturally led to a strong fractionation, in which more than half of the planet's heat producing element content may be present in the crust²⁸. In this so-called stagnant lid regime, the crust further acts as a physical boundary that thermally insulates the interior, dictating the rate at which Mars cools in time⁷,¹⁸. Thus, the structure and properties of the crust are key to decipher the geodynamic evolution of Mars throughout its geologic history.

We simulate Mars's three-dimensional (3D) thermal evolution and interior viscous flow using the geodynamic code GAIA²⁹ and including the latest constraints from InSight on the interior structure, namely core size and crustal thickness²⁵,³⁰ (Methods). Extrapolating InSight local crustal thickness measurements to planetary scale using gravity and topography data is non-unique, and different models accounting for variations in crustal density are tested²⁶. Crustal heat production is assumed to be laterally constant, but is ranged from the Gamma Ray Spectrometer measured surface average of 49 pW kg⁻¹ (ref. 28), to up to twice that value. In total, 84 global geodynamic models are created, among which a subset of models covering a wide range of interior viscosity and temperature are retained for further analysis (Methods, Extended Data Table 1 and Extended Data Fig. 2).

The present-day interior structure beneath the north polar cap is represented for each model by averaging all cells poleward of 60° N in the geodynamic models. In all simulations, the interior shows a stratified viscosity structure with a thick stagnant lid (thicker than 350 km) on top of a lower viscosity mantle (10²¹–10²⁴ Pa s) with upper mantle (100–400 km) thermal gradients of 2–3 K km⁻¹, consistent with earlier work¹⁶,³¹,³². For reference, our subset of interior structures are plotted alongside pre-InSight models¹⁸ (Fig. 2). Models with a lower or higher viscosity would require crustal thicknesses thinner or larger than constrained by InSight²⁶, respectively. In all models, a larger crustal thickness and a higher concentration of heat producing elements in the crust result in a cooler and stiffer mantle. In absence of additional constraints, all selected models are equally possible.

## Ice history and viscoelastic relaxation

Global climate models considering Mars's recent orbital history³³ suggest that the current north polar cap formed over the past 4 million years, with ice-sheet thickness increasing in response to a gradual decline in polar insolation²³. A geologically recent formation is also supported by crater counting statistics and stratigraphic relationships that indicate ages of only a few million years, albeit with significant uncertainties²². Although these studies suggest a young age for the current north polar cap, its precise age and formation history remain uncertain¹³. Former ice caps may have also existed²³,³⁴, but the planet's orbital and climate history cannot be uniquely determined beyond the last several tens of millions of years³³, making it difficult to infer longer term ice histories. If the ice loading time is shorter than or comparable to the adjustment time of the lithospheric mantle (tens to hundreds of millions of years), the equilibrium may not be reached resulting in an ongoing viscoelastic downward deflection. This effect would prominently affect all estimates of the strength of the elastic lithosphere and planetary thermal state, in which viscous relaxation is not addressed adequately⁵,⁶,³⁵.

Viscoelastic relaxation is investigated by computing the time-variable and wavelength-dependent load Love number *h*′ using the ALMA code, which solves the momentum equation of a spherically symmetric stratified viscoelastic planet⁸ (Methods). Our calculations use a Maxwell rheology, as it is established to describe long-term glacial isostatic adjustment processes¹,². Because the age and long-term loading history of the north polar cap are poorly known, we assume a linear increase in thickness over time, leading to the current state, with tested final ages ranging from 100 thousand years ago (ka) to 1 billion years ago (Ga). Although nonlinearities in the ice-accumulation history are expected³⁶, these would have a negligible effect on present-day deformations (Extended Data Fig. 3). The effect of former north polar caps on the current polar deformations is also explored (Methods). Love numbers are computed for each of our interior models up to spherical harmonic degree 50 and at 300 times spaced evenly on a log scale from 100 ka to 1 Ga. Changes in the seasonal cap whose load is three orders
<figure><img src="image_4.png" /><figcaption>**Fig. 2** | Present-day average interior structure of the northern regions (poleward of 60° N). a, b, Coloured lines show the viscosity (a) and temperature (b) structure of the interior models that match InSight constraints. Line colour</figcaption></figure>

reflects the average mantle temperature, which is related to the thickness of the crust and its enrichment in heat producing elements (HPE). Detailed model description is shown in Extended Data Fig. 1 and Extended Data Table 1.

of magnitude lower than the perennial cap³⁷ and that are insensitive to the interior viscosity⁴ are not considered.

## Constraints from radar observations

The present-day north polar cap has a thickness of about 3.1 km and bulk density of $1,200 \pm 300 \text{ kg m}^{-3}$, as constrained from radar analyses⁶,³⁴, and together with our specified ice loading history these are used to predict the time-variable polar deformations (Methods). For each model, the time-integrated deformations are compared to radar observations of the undeformed northern cap basement from Mars Advanced Radar for Subsurface and Ionosphere Sounding (MARSIS¹⁴) and SHAIIow RADar (SHARAD⁵). The misfit is evaluated using a root-mean-square (r.m.s.) and our viscoelastic and radar analyses consider realistic combinations of permittivity and density for a mixture of dust, water and dry ice⁶ (Methods). Accounting for a wide range of uncertainties affecting radar measurements, including instrument range resolution, surface roughness and the uncertainty in the estimation of the preloading polar cap basement, the r.m.s. misfit should be no more than 266 m for MARSIS and 175 m for SHARAD (Methods). For all models, the flexure beneath the polar cap increases as a function of time (or polar cap age) and is associated with an increase in the cumulative radar misfit (sum of SHARAD and MARSIS misfits, Fig. 3a). At all investigated time steps, flexure is greater for lower viscosity interiors due to these models reaching equilibrium more rapidly.

## Mars's time-variable gravity field

Ongoing deformations affect the planetary gravity field³⁸ and in particular the long-wavelength zonal degrees 2 and 3. For each model, we compute the present-day yearly surface deformation rate linked to mantle flow. In all cases, the younger the polar cap, the higher the present-day yearly deformation rate, which is related to the time-variable viscoelastic response of the interior (Fig. 3b). These values are compared to observational constraints from orbital tracking of the Mars Global Surveyor, Mars Odyssey and Mars Reconnaissance Orbiter missions⁹, which show that the planetary degree 2 and 3 zonal gravity field coefficients have marginally increased over 8 Martian years of spacecraft tracking (Extended Data Fig. 4). This variation indicates that the increase in the gravitational potential associated with long-term ice accumulation is higher than the decrease in gravitational potential from downward deflection. Based on climate

models, ice-accumulation rates have been constrained to range from 0.5 to 1.1 mm per year at the north pole²³,³⁹. Considering density ratios between the accumulating ice and flowing interior mantle, the above minimum accumulation rate implies that ongoing downward deformation in the northern regions should be less than 0.13 mm per year to match the observed time-variable gravity field (Methods).

## Low polar strain rates

If the northern strain rates were too large, crustal failure would occur, potentially leading to the detection of a marsquake by InSight¹⁰. Therefore, the polar deformations can further be constrained by comparing the geodynamic and seismic moment rates⁴⁰, as has been done on Earth⁴¹. Based on predicted magnitudes from the InSight catalogue, a 3.8 to 3.9 magnitude event originating from the polar regions (75–90° N) could be detectable¹⁰. Thus, the non-observation of such marsquakes during the InSight mission can be used to invert for the maximum northern lithosphere strain rate (Methods). This approach shows that the viscoelastic strain rate should be less than $1.84 \times 10^{-18} \text{ s}^{-1}$ to avoid inducing a 3.8 magnitude marsquake poleward of 75° N.

Together, the above constraints are used to limit our multidimensional parameter space encompassing the interior structure, age and bulk properties of the north polar cap. A model is deemed acceptable if the r.m.s. misfit between the present-day deformation and radar measurements is less than both 266 and 175 m for MARSIS and SHARAD, respectively, if the combination of tested density and dielectric constant produce an existing mixture of ice and dust, and if both the deformation and strain rates do not exceed 0.13 mm per year and $1.84 \times 10^{-18} \text{ s}^{-1}$, respectively. If the polar cap were older, it would have had more time to adjust leading to a vertical displacement below the ice cap larger than observed, whereas if it were younger, deformation rates would be too large and inconsistent with InSight seismic moment rate and Mars's gravity field (Fig. 3a,b). In our inversion, radar measurements provide the tightest constraints on the interior viscosity and limit the maximum age of the polar cap, whereas the time-variable gravity provides a lower bound on the polar cap age. The InSight moment rate establishes a weaker constraint on the polar cap minimum age, but is not in contrast with our gravity analyses.

## Highly viscous mantle and young ice cap

Only three members of our geodynamically based model ensemble were found to be consistent with the small amount of deformation
<figure><img src="image_5.png" /><figcaption>**Fig. 3** | Observational constraints on glacial isostatic adjustment from radar measurements, time-variable gravity and InSight moment rate. a, North polar cap age versus cumulative misfit of viscoelastic and measured deformation from SHARAD and MARSIS data for all models. b, Time-variable deformation and strain rate for all models. c, Polar cap age and maximum</figcaption></figure>

deformation beneath the north polar cap for accepted models. The colours are the same as in Fig. 2 and the black square in **b** and **c** shows values for the 4 million years old polar cap loading model in Fig. 4. Further information on models is given in Extended Data Fig. 1 and Extended Data Table 1. Max., maximum. Ma, million years ago.

constrained by radar data, time-variable gravity field and InSight seismicity (Fig. 3c). These models are characterized by a large present-day mantle viscosity with volumetric average values of about $2-6 \times 10^{22}$ Pa s from 500 km down to the core and thick average crust (thicker than 40 km, Figs. 3 and 4). Models with a lower viscosity would predict too much north polar deformation, whereas models with a higher viscosity would require crustal thicknesses larger than measured by InSight²⁶. Although our inversion constrains the interior structure of the northern regions, it is comparable to the global average, albeit marginally colder (Extended Data Fig. 5). In our accepted models, the mantle is highly depleted in radiogenic elements compared to the crust, the latter of which must have a heat production rate of at least 61 pW kg⁻¹ for our 69-km average crustal thickness model, and up to 98 pW kg⁻¹ for our 43-km crust model. These values are 1.25 and up to two times greater than the surface average²⁸ and indicates that more than 90% of the planet's heat producing element content is present in the crust. This shows that the lower crust has a heat producing element concentration and composition different from that of the surface, as previously suggested³⁵,⁴²,⁴³. In agreement with previous work²⁶,⁴³, our estimated average crustal thickness implies that the crust-mantle interface at the InSight location corresponds to the deepest observed reflector (that is, the three-layered model³⁰). Our low preferred mantle temperature further agrees with that inferred from the postolivine transition detected at 1,000 km by InSight⁴⁴.

Our models limit the age of the north polar cap load to be 1.7–12.0 Ma, which is consistent with independent global climate models suggesting ages of a few million years²³. If the age of the polar cap or interior structure of Mars were known, our framework would enable to place precise constraints on either one, potentially allowing to refine the long-term orbital history of the planet and ice-accumulation history of the north polar cap. Based on our allowed deformation rate and the

non-observation of high-magnitude marsquakes poleward of 75° N, the ratio between the northern lithosphere deformation rate and seismic moment rate (that is, the seismic efficiency) is found to be at most 0.3. This value is more than two times lower than found using InSight data alone¹⁰. Such a low seismic efficiency helps explain why InSight has detected fewer marsquakes than expected²⁴ and may be related to large crustal porosity or the presence of volatiles in the crust.

<figure><img src="image_6.png" /><figcaption>**Fig. 4** | Constrained present-day average interior structure of the northern regions (poleward of 60° N) and polar deformations. **a**, Accepted interior viscosity and temperature ranges of the northern regions. **b**, Example of the present-day deformation and subsidence rate for a polar cap age of 4 million years for our 47–3.0_1.25λ model (Extended Data Table 1).</figcaption></figure>
Given the long interior response to loading, we have considered models accounting for the potential existence of former north polar caps. Simulations assuming ice cap formation following periodic obliquity cycles only moderately increased the north polar deformation as a result of continuing adjustment of the ancient uncompensated loads (Extended Data Figs. 6–8). This would not affect our inferred mantle viscosities or heat productions, but would lower the maximum current north polar cap age by a few million years. One way to help reduce the ongoing deformation, and thus our inferred interior viscosity, would require the present-day polar cap to have formed on the sedimentary-infilled flexural trough of a former ice cap. In that framework, the ongoing downward deformation would compete with the rebound related to the past ice cap. However, although such sedimentary infilling may be related to the basal unit beneath the central portion of the polar layered deposits³⁴, it is not seen beneath the Gemina Lingula region near 80° N (Fig. 1). Therefore, this competing effect alone cannot account for the negligible flexure observed in that region.

Thus far, analyses of glacial isostatic adjustment have been limited to Earth applications because of a lack of observational data. Our work shows that there is opportunity to study this process on another planet and to unveil the interior structure, geodynamic history and long-term orbital evolution of Mars. Furthermore, our analysis sheds light on viscosity variations in the interior of a stagnant lid body, for which long-term geodynamic constraints are scarce. The ongoing subsidence rates of 10⁻² to 10⁻¹ mm per year predicted by our models affect the time-variable planetary gravity field, a signal that could be directly measured by future space missions to Mars¹¹,¹² thereby bridging the gap in geodynamic modelling and interior exploration between Mars and Earth.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-024-08565-9.

1. Peltier, W. R. The impulse response of a Maxwell Earth. *Rev. Geophys.* **12**, 649–669 (1974).

2. Sabadini, R., Vermeersen, B. & Cambiotti, G. *Global Dynamics of the Earth* (Springer, 2016).

3. Greve, R., Klemann, V. & Wolf, D. Ice flow and isostasy of the north polar cap of Mars. *Planet. Space Sci.* **51**, 193–204 (2003).

4. Wagner, N. L., James, P. B., Ermakov, A. I. & Sori, M. M. Evaluating the use of seasonal surface displacements and time-variable gravity to constrain the interior of Mars. *J. Geophys. Res. Planets* **129**, e2023JE008053 (2024).

5. Phillips, R. J. et al. Mars North Polar Deposits: stratigraphy, age, and geodynamical response. *Science* **320**, 1182–1185 (2008).

6. Broquet, A., Wieczorek, M. A. & Fa, W. Flexure of the lithosphere beneath the north polar cap of Mars: implications for ice composition and heat flow. *Geophys. Res. Lett.* https://doi.org/10.1029/2019GL086746 (2020).

7. Plesa, A.-C. et al. in *Geophysical Exploration of the Solar System* Vol. 63 of series *Advances in Geophysics* (eds Schmelzbach, C. & Stähler, S. C.) 179–230 (Elsevier, 2022).

8. Melini, D., Saliby, C. & Spada, G. On computing viscoelastic Love numbers for general planetary models: the ALMA3 code. *Geophys. J. Int.* **231**, 1502–1517 (2022).

9. Genova, A. et al. Seasonal and static gravity field of Mars from MGS, Mars Odyssey and MRO radio science. *Icarus* **272**, 228–245 (2016).

10. Knapmeyer, M. et al. The global seismic moment rate of Mars after event S1222a. *Geophys. Res. Lett.* **50**, e2022GL02296 (2023).

11. Wörner, L. et al. MaQuls—concept for a Mars quantum gravity mission. *Planet. Space Sci.* **239**, 105800 (2023).

12. Genova, A. ORACLE: a mission concept to study Mars' climate, surface and interior. *Acta Astronaut.* **166**, 317–329 (2020).

13. Byrne, S. The polar deposits of Mars. *Annu. Rev. Earth Planet. Sci.* **37**, 535–560 (2009).

14. Picardi, G. et al. Radar soundings of the subsurface of Mars. *Science* **310**, 1925–1928 (2005).

15. Khan, A. et al. Upper mantle structure of Mars from InSight seismic data. *Science* **373**, 434–438 (2021).

16. Pou, L. et al. Tidal constraints on the Martian interior. *J. Geophys. Res. Planets* **127**, e2022JE007291 (2022).

17. Kiefer, W. S. & Li, Q. Mantle convection controls the observed lateral variations in lithospheric thickness on present-day Mars. *Geophys. Res. Lett.* https://doi.org/10.1029/2009GL039827 (2009).

18. Plesa, A.-C. et al. The thermal state and interior structure of Mars. *Geophys. Res. Lett.* **45**, 12,198–12,209 (2018).

19. Hauber, E., Brož, P., Jagert, F., Jodlowski, P. & Platz, T. Very recent and wide-spread basaltic volcanism on Mars. *Geophys. Res. Lett.* https://doi.org/10.1029/2011GL047310 (2011).

20. Voigt, J. R. C. et al. Revealing Elysium Planitia's young geologic history: Constraints on lava emplacement, areas, and volumes. *J. Geophys. Res. Planets* **128**, e2023JE007947 (2023).

21. Broquet, A. & Andrews-Hanna, J. Geophysical evidence for an active mantle plume underneath Elysium Planitia on Mars. *Nat. Astron.* **7**, 160–169 (2023).

22. Tanaka, K. Geology and insolation-driven climatic history of Amazonian north polar materials on Mars. *Nature* **437**, 991–994 (2005).

23. Levrard, B., Forget, F., Montmessin, F. & Laskar, J. Recent formation and evolution of northern Martian polar layered deposits as inferred from a global climate model. *J. Geophys. Res. Planets* https://doi.org/10.1016/10.1029/2006JE002772 (2007).

24. Banerdt, W. et al. Initial results from the InSight mission on Mars. *Nat. Geosci.* **13**, 183–189 (2020).

25. Stähler, S. C. et al. Seismic detection of the Martian core. *Science* **373**, 443–448 (2021).

26. Wieczorek, M. A. et al. InSight constraints on the global character of the Martian crust. *J. Geophys. Res. Planets* **127**, e2022JE007298 (2022).

27. Taylor, G. J. The bulk composition of Mars. *Geochemistry* **73**, 401–420 (2013).

28. Hahn, B. C., McLennan, S. M. & Klein, E. C. Martian surface heat production and crustal heat flow from Mars Odyssey Gamma-Ray Spectrometry. *Geophys. Res. Lett.* https://doi.org/10.1029/2011GL047435 (2011).

29. Hüttig, C., Tosi, N. & Moore, W. B. An improved formulation of the incompressible Navier-Stokes equations with variable viscosity. *Phys. Earth Planet. Inter.* **220**, 11–18 (2013).

30. Knapmeyer-Endrun, B. et al. Thickness and structure of the martian crust from InSight seismic data. *Science* **373**, 438–443 (2021).

31. Samuel, H., Lognonné, P., Panning, M. & Lainey, V. The rheology and thermal history of Mars revealed by the orbital evolution of Phobos. *Nature* **569**, 523–527 (2019).

32. Khan, A., Huang, D. & Durán, C. E. A. Evidence for a liquid silicate layer atop the Martian core. *Nature* **622**, 718–723 (2023).

33. Laskar, J. et al. Long term evolution and chaotic diffusion of the insolation quantities of Mars. *Icarus* **170**, 343–364 (2004).

34. Nerozzi, S., Ortiz, M. R. & Holt, J. W. The north polar basal unit of Mars: an Amazonian record of surface processes and climate events. *Icarus* **373**, 114716 (2022).

35. Ojha, L., Karimi, S., Lewis, K. W., Smrekar, S. E. & Siegler, M. Depletion of heat producing elements in the Martian mantle. *Geophys. Res. Lett.* **46**, 12756–12763 (2019).

36. Smith, I. B. et al. An ice age recorded in the polar deposits of Mars. *Science* **352**, 1075–1078 (2016).

37. Xiao, H. et al. Spatio-temporal level variations of the Martian seasonal north polar cap from co-registration of MOLA profiles. *J. Geophys. Res. Planets* **127**, e2021JE007158 (2022).

38. Tamisiea, M. E., Mitrovica, J. X. & Davis, J. L. GRACE gravity data constrain ancient ice geometries and continental dynamics over Laurentia. *Science* **316**, 881–883 (2007).

39. Laskar, J., Levrard, B. & Mustard, J. Orbital forcing of the Martian Polar Layered Deposits. *Nature* **419**, 375–377 (2002).

40. Phillips, R. J. Appendix B. Expected rates of marsquakes. In *Scientific Rationale and Requirements for a Global Seismic Network on Mars* LPI Technical Report 91-02 (NASA, 1991).

41. Ward, S. N. On the consistency of earthquake moment rates, geological fault data, and space geodetic strain: the United States. *Geophys. J. Int.* **134**, 172–186 (1998).

42. Thiriet, M., Michaut, C., Breuer, D. & Plesa, A.-C. Hemispheric dichotomy in lithosphere thickness on Mars caused by differences in crustal structure and composition. *J. Geophys. Res. Planets* **123**, 823–848 (2018).

43. Kim, D. et al. Global crustal thickness revealed by surface waves orbiting Mars. *Geophys. Res. Lett.* **50**, e2023GL013482 (2023).

44. Huang, Q. et al. Seismic detection of a deep mantle discontinuity within Mars by InSight. *Proc. Natl Acad. Sci. USA* **119**, e2204474119 (2022).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_7.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2025
## Constraints on the deformation beneath the north polar cap

The deformation beneath the north polar cap is obtained by comparing estimates on the polar cap thickness from elevation and radar data following ref. 6. In that work, the polar cap thickness was estimated by subtracting the observed ice cap elevation from a preloading surface constructed using elevation data far from the polar cap and unaffected by flexure (Extended Data Fig. 9). MARSIS radar analyses were performed at 213 locations spatially scattered across the polar deposits. All available MARSIS radargrams close to each location were investigated, and the reflections arising from the icy surface and the ice-substratum interface were visually identified. Using the same framework, analyses of SHARAD radargrams were performed at those 213 locations. Because of its higher frequency, SHARAD does not penetrate through the sand-rich basal unit and we here only retain locations where the ice-crust interface is observed (n = 78).

The radar thickness depends on the dielectric properties of the icy materials, which are not uniquely known for the entire north polar cap. On the other hand, viscoelastic deformations depend on both the interior structure and density of the loading materials. Therefore, we compare our simulated viscoelastic deformations to radar measurements using an r.m.s. misfit function for all interior models considering a real dielectric constant range of 2.5 to 3.5 and an ice density of 920 to 1,500 kg m⁻³ covering the possible presence of water and dry ice and dust in the ice cap⁵,⁶. We limit our ice density and dielectric constant parameter space by considering possible mixtures of dust, water and dry ice using a three-component Maxwell-Garnett mixing model⁶. The density of water and dry ice is considered to be 920 and 1,560 kg m⁻³, respectively, and the dust density is assumed to range from 2,200 to 3,400 kg m⁻³. The dielectric constants of these same materials are 3.0, 2.5 and 6.0, respectively⁵,⁶. Because of these, a model with an ice density of 920 kg m⁻³ (pure water ice), but with a dielectric constant of 3.5, can be safely ruled out as it does not correspond to any existing mixture of ice and dust in the north polar cap.

By accounting for uncertainties linked to the range resolution of MARSIS, surface roughness at the scale of the Fresnel zone of MARSIS and the uncertainty when estimating the preloading surface, an upper limit on the r.m.s. was estimated to be 266 m (ref. 6). Using the same approach and because of the higher frequency and resolution of SHARAD, a lower maximum allowed misfit of 175 m is given to the SHARAD comparison. Thus, a model is deemed acceptable if the r.m.s. misfit for MARSIS and SHARAD is no more than 266 and 175 m, respectively.

For context, Fig. 1b shows a flexure map beneath the north polar cap that was obtained by interpolating the radar-derived flexure using a 400-km moving window to obtain a smooth basement map and to get rid of short-wavelength uncertainty in elevation and radar-derived polar cap thickness. For comparison, a similar map using a dielectric constant of 2.75 instead of 3.0 is shown in Extended Data Fig. 9. In both cases, the present-day deformation beneath the north polar cap is no more than a few hundreds of metres.

## Thermal evolution modelling

We have run 84 Martian 3D geodynamic models considering a wide range of parameters using the GAIA code⁷,¹⁸,²⁹. Therein, conservation equations of mass, linear momentum and thermal energy, are solved from 4.5 GPa to the present day given a set of model parameters. Key parameters that control the present-day thermal state of the interior are core radius, crustal thickness and radiogenic heat production, some of which have been recently constrained by InSight²⁵,³⁰. Our geodynamic simulations use the 3D structure of the crust, as derived from gravity, topography and InSight data²⁶. Absolute viscosities can be obtained using the Arrhenius viscosity law and considering reference values for the viscosity, pressure and temperature⁷.

All tested models use the following naming convention: $d_c^{\text{InSight}} - \rho_{\text{south}}(-\rho_{\text{north}})$, where $d_c^{\text{InSight}}$ is the crustal thickness at the InSight landing site in kilometres and $\rho_{\text{north}}$ and $\rho_{\text{south}}$ are the bulk density of the northern and southern hemisphere crust in grams per cubic centimetre (Extended Data Fig. 1 and Extended Data Table 1). If only one crustal density, $\rho$, is indicated, this value is assumed to be constant across the planet. If added, $\lambda$ provides the crustal heat producing element enrichment factor with respect to the nominal Gamma Ray Spectrometer measured average (below). Models with additional annotations have higher activation energy ($E = 325 \text{ kJ mol}^{-1}$), initial temperature ($T_{\text{init}} = 1,850 \text{ K}$), temperature difference across the mantle ($\Delta T = 2,205 \text{ K}$) and have a $\eta_{\text{jump}}$ 25-fold mantle viscosity jump compared to the reference values used in most other cases (that is, $E = 300 \text{ kJ mol}^{-1}$, $T_{\text{init}} = 1,650 \text{ K}$, $\Delta T = 2,000 \text{ K}$, no viscosity jump, see Extended Data Table 1 for further details). The effect of these parameters within their allowed possible ranges is minor compared to crustal thickness (Fig. 2).

In our nominal models, crustal heat production is laterally constant and equal to the surface average of 49 pW kg⁻¹, estimated from Gamma Ray Spectrometer data²⁸. However, because the spectrometer instrument is only sensitive to the first upper tens of centimetres of the surface⁴⁵, different composition and heat production may exist at greater depth. Therefore, we also consider models in which the bulk crust has a heat production different from the near surface, with an enrichment factor, $\lambda$, of 0.83 and up to 2. The upper bound leads to heat production of 98 pW kg⁻¹, which is slightly higher than the highest measured crustal heat production of 75 pW kg⁻¹ (ref. 28), making it a reasonable upper range. Because the number of radiogenic elements in a planetary interior is finite, some models with a thick crust cannot reach crustal enrichment factors of 2, as the mantle is already fully depleted. Thus, our approach establishes an upper limit for Mars's mantle viscosity (and the lowest possible temperature) within our geodynamic evolution framework. Whereas higher viscosities may be obtained using a higher reference viscosity, these models would be mostly conductive and be unable to explain ongoing volcanism on Mars.

Geologic observations and geophysical models indicate the presence of melts in the Martian interior¹⁹⁻²¹. All of our models are able to predict melting at the present day when considering the present-day Martian solidus as estimated in ref. 46. From our initial set of 84 geodynamic models, we select a subset of 27 end-member models with a wide range of viscosity structure for our viscoelastic analyses (Extended Data Table 1).

## Thermal profile to lithosphere elastic thickness

A thermal profile can be converted to an elastic lithosphere by setting the bending moment of an elastic plate equal to that of the bending stresses in a more realistic rheology that considers fracturing and viscous flow⁴⁷ as

$$ \frac{E_y T_c^2 K}{12(1 - v^2)} = \int_0^{T_m} \sigma_{\text{YSE}}(z)(z - z_n) dz, \quad (1) $$

where the left term corresponds to the analytic integration of the bending moment of an elastic plate with thickness $T_c$, and in which $v$ is Poisson's ratio, $E_y$ is Young's modulus and $K$ is the lithosphere curvature. The right term computes the bending moment limited by the yielding strength of the lithosphere ($\sigma_{\text{YSE}}$), and where $z$ is the depth, $z_n$ the depth of the neutral plane and $T_m$ is the mechanical thickness⁴⁷. The yield strength of the lithosphere strongly depends on temperature variations and includes faulting and viscous stresses. The mechanical thickness is defined as the depth where bounding stresses are close to zero, here approximated to 50 MPa (refs. 5,47).

Thermal profiles from our interior models are converted to an equivalent elastic thickness using the TeHF code⁴⁸ that solves the above equation. We consider a wide range of model parameters including, a dry and wet diabase rheology for the crust, a dry and wet olivine rheology
for the mantle. For each rheology, we consider strain rates of $10^{-17}$ to $10^{-20}$ s$^{-1}$ consistent with the timescale of mantle convection and polar cap formation¹⁸ (Fig. 3). The curvature of the lithosphere is also varied from $10^{-9}$ to $10^{-11}$ m$^{-1}$ (ref. 6). For all models, this conversion approach leads to elastic thicknesses less than 270 km, which is inconsistent with constraints from elastic flexure indicating elastic thicknesses more than 330 km (refs. 5, 6, Extended Data Fig. 1). This demonstrates that previous elastic modelling must be revisited to consider the interior transient viscoelastic response.

## Calculation of viscoelastic deformation

The vertical displacement at time $t_{end}$, colatitude and longitude ($\theta, \phi$), are determined by the convolution of the spectral representation of the load potential change over the considered time interval with the transfer function represented by the load Love number $h'$ as

$$W(t_{end}, \theta, \phi) = \int_{0}^{t_{end}} \sum_{l=0}^{l_{max.}} \sum_{m=-l}^{l} \frac{V_{lm}(t)Y_{lm}(\theta, \phi)}{g_{0}} h'_{l}(t_{end} - t)dt, \quad (2)$$

where $g_0$ is the gravitational acceleration at the surface, $l$ and $m$ are the degree and angular order, $Y_{lm}$ are real spherical harmonics and $l_{max.}$ is the maximum spherical harmonic degree used for calculations. The change in the load potential of the north polar cap ($V_{lm}$) is calculated assuming a mean bulk density ranging from 920 to 1,500 kg m$^{-3}$ and using the polar cap thickness derived from elevation data and radar measurements$^{5,6}$. For simplicity and because of a lack of observational data, we solely increase the polar cap thickness over time, based on the current morphology of the north polar cap. We do not consider any radial expansion or growth.

Love numbers are computed based on the internal structures predicted by thermal models (Extended Data Fig. 2) using a Maxwell rheology and up to $l_{max.} = 50$, which is largely sufficient to resolve the polar cap loading (roughly degree 8, wavelength of 1,250 km). We show that the effect of using a transient rheology such as Andrade is minor compared to the uncertainty in mantle viscosity structure, and note that it does not change the time behaviour for the considered long-term process (Extended Data Fig. 10). The ALMA code assumes an incompressible interior. Compressibility affects Love number calculation, leading to slightly higher deformations. This effect is considered as second order when compared to the uncertainty in Mars's interior structure.

Each interior model from our thermal evolution simulation is discretized into 68 constant-thickness layers with a given density, viscosity and rigidity, from the core to the surface. This number of layers is chosen to optimize the calculation of the load Love number using ALMA. The interior rigidity is computed using Perple_X (ref. 49) together with the interior temperature, pressure and TAYAK mineralogy$^{7,26,27}$.

This study uses average one-dimensional interior profiles to constrain viscoelastic deformations beneath the north polar cap. Although future work may achieve larger accuracy by modelling 3D deformations from a 3D interior structure and polar cap load, several mitigating factors should be noted. Our 3D geodynamic models indicate that Mars's northern regions show homogeneous properties (rigidity, viscosity, temperature) comparable to the global average (Extended Data Fig. 5). This indicates that lateral variations in the interior structure have minor effects on the estimated deformations. In particular, we expect these effects to be markedly lower than the current uncertainty in Mars's interior structure (Fig. 2). Furthermore, the axisymmetric shape of the north polar cap load reduces the influence of the polar cap geometry on deformation estimates. For these reasons, differences between a full 3D viscoelastic loading model and our current model are expected to be small.

## Polar cap loading history and past ice caps

On Mars, polar caps have grown and decayed following the planet's obliquity cycles$^{23}$. It is not uniquely known how obliquity has varied

in the planet's recent history (500 million years ago), as orbital models become chaotic as one goes backward in time$^{33}$. Yet, it is generally thought that past polar caps existed in the geologically recent history$^{23,34}$, and these may affect the present-day deformation in the northern regions.

To model this process, we have constructed a loading history in which north polar caps grow and decay. Our first model starts at 430 Ma with no polar cap, and builds up a 3.2-km-thick polar cap in 5 Myr following a half cycle of a cosine function, removes that polar cap in 5 Myr and then repeats with a new cycle. The radius of the edge of the polar cap is constant. The interior response associated with the ice cap periodic growth and decay is modelled over the full process duration using the viscoelastic formalism described above. The two end-member viscosity structures from our ensemble of geodynamic models are tested with this time-loading history. In both cases, the flexure values at the present day are 35 and 45 m larger than when neglecting this past history for the low- and high-viscosity models, respectively (Extended Data Fig. 6).

We have also tested these same interior models against the ice-accumulation history of ref. 23. This model starts at 10 Ma with a past north polar cap that disappears at about 8 Ma. After a time interval with no north polar cap, a new polar cap forms with a near linear increase in thickness from 4 Ma to the present day (Extended Data Fig. 7). A key unknown is what was the state of the 10 Ma ice cap and whether it was long-standing and at equilibrium or more recently formed. In one model, we assume that the previous ice cap was present over the last 500 Ma and thus at equilibrium (Extended Data Fig. 7) and in a second model, we assume that it formed at 14 Ma (Extended Data Fig. 8). In the case where the ice cap was long lasting, the present-day deformation increased by up to 70 m, whereas it only increased by a few metres for a young former ice cap.

The above tests show that the deformation beneath the north pole can only be increased when considering the effect of ancient ice caps. Whereas this effect would reduce our inferred maximum age for the polar cap by a few million years ago, it would not affect our constraints on Mars's interior structure.

## Constraints on polar deformation from time-variable gravity

Previous work inverted for Mars's static and time-varying gravity field on the basis of tracking of Mars Global Surveyor, Mars Odyssey and Mars Reconnaissance Orbiter missions over 8 Martian years$^9$. Here we build on that work by analysing residual trends in $C_{20}$ and $C_{30}$ after the zonal harmonics are corrected for seasonal variability. In both cases, the sign of these coefficients is negative, indicating the planet's gravitational oblateness and north-south asymmetry. We correct the time-varying signal of these zonal harmonics for annual, semi-annual and tri-annual variations, as well as half (5.5 years) and full solar cycle periods (11 years, ref. 9). The function used to fit the both $C_{20}$ and $C_{30}$ time-variable coefficients is in the form of

$$f(t) = P_1 \sin\left(\frac{2\pi}{T_1} t\right) + P_2 \cos\left(\frac{2\pi}{T_1} t\right) + P_3 \sin\left(\frac{2\pi}{T_2} t\right) \\ + P_4 \cos\left(\frac{2\pi}{T_2} t\right) + P_5 \sin\left(\frac{2\pi}{T_3} t\right) + P_6 \cos\left(\frac{2\pi}{T_3} t\right) \\ + P_7 \sin\left(\frac{2\pi}{T_4} t\right) + P_8 \cos\left(\frac{2\pi}{T_4} t\right) + P_9 \sin\left(\frac{2\pi}{T_5} t\right) \\ + P_{10} \cos\left(\frac{2\pi}{T_5} t\right), \quad (3)$$

where $P_1$ to $P_{10}$ are the best-fit periodic coefficients, $t$ is time in seconds past J2000, $T_1 = 365 \times 86,400 \times 1.880894$, $T_2 = 365 \times 86,400 \times 0.940447$, $T_3 = 365 \times 86,400 \times 0.626965$, $T_4 = 11 \times 86,400 \times 365$ and $T_5 = 11 \times 86,400 \times 365/2$.

After this correction, the residuals show trends with slopes of $1.5 \pm 1.6 \times 10^{-18}$ and $2.5 \pm 3.3 \times 10^{-19}$ s$^{-1}$ for the $C_{20}$ and $C_{30}$ coefficients,
respectively (Extended Data Fig. 4 and Extended Data Table 2). Given the negative sign of C₂₀ and C₃₀, positive trends suggest polar ice accumulation and that the south polar cap is accumulating less ice compared to the north polar cap. For those time-variable coefficients the Pearson correlation coefficients (r values) are 0.13 and 0.06, respectively. We further test the null hypothesis that the residuals have no correlation with time using Wald's test. For the C₂₀ and C₃₀ coefficient residuals, we obtain P values of 0.35 and 0.46, respectively, which are both higher than typical significance levels of 0.05. Together, these values indicate a positive slope, although not statistically significant.

After 1 Martian year, our analysis indicates a ΔC₂₀ value of 4.7 × 10⁻¹¹ and ΔC₃₀ of 7.7 × 10⁻¹². Together, these can be used to derive the north polar gravitational acceleration as

$$a_m \approx 3 \frac{GM}{R^2} \Delta C_{20} + 4 \frac{GM}{R^2} \Delta C_{30}, \quad (4)$$

where G is the gravitational constant, M and R are Mars's mass and radius. Based on this equation, our measured ΔC₂₀ and ΔC₃₀ imply north polar gravitational acceleration of 6.4 × 10⁻⁵ mGal. Using the Bouguer plate formula, 2πGph, and a water ice density of 920 kg m⁻³, the above gravity acceleration suggests a maximum ice-accumulation rate of 1 mm per year when not accounting for isostatic adjustment. This value is consistent with predicted rates of 0.50 to 1.06 mm per year by climate models²³,³⁹, indicating that gravity residuals provide information on north polar processes.

Our models predict that glacial isostatic adjustment is ongoing, with downward deformation rates of 10⁻⁴ to 4 mm per year (Fig. 3). In our framework, glacial isostatic adjustment pushes the mantle downward such that the gravitational signature of this process should be scaled by mantle density (3,500 kg m⁻³) and has an effect opposite to ongoing ice accumulation that scales with water ice density (920 kg m⁻³). Given the overall purity of the north polar layered deposits⁵, dust accumulation is not accounted for here. Erosion rates that have been measured to be small, 2 × 10⁻⁴ mm per year (ref. 50), are also neglected. Owing to its higher elevation, the south polar cap is expected to only accumulate small amount of dry ice⁵¹, which is here neglected.

Because of a lack of a clear trend in our residual gravity analyses, we here assume that the gravitational potential from ongoing subsidence cannot be higher than from the polar cap accumulation rate. This indicates that ongoing subsidence rates must be at most around 26% (920/3,500) that of ice accumulation. Ice-accumulation estimates from previous work therefore limit the subsidence rate to less than 0.13 mm per year (ref. 39) or less than 0.28 mm per year (ref. 23) and we here use the minimum of these two values to limit our parameter space (Fig. 3). Choosing the other value has no effect on our derived viscosity structure, but would provide a lower limit to the polar cap age of 0.9 Ma instead of 1.7 Ma.

## Moment rate inversion from InSight

The InSight seismometer has detected dozens of marsquakes over its 4 years of activity²⁴, but none from the north polar regions (75–90° N). Analyses of InSight data have shown that a 3.8 to 3.9 magnitude event originating from these northern regions could have been detected¹⁰. Because the seismic moment rate depends on the strain rate experienced by the seismically competent lithosphere⁴⁰,⁴¹, it is possible to invert for a maximum strain rate in the northern regions based on this non-observation. The maximum strain rate is obtained by comparing the seismic moment rate (Ṁ_seismic), as measured by InSight, to a geodetic moment rate (Ṁ_geodetic) obtained from the Kostrov equation⁴⁰.

Mars's seismic moment rate based on the InSight catalogue can be defined in terms of the b value¹⁰ as

$$\dot{M}_{\text{seismic}} = \frac{\Gamma(2 - 2b/3)}{n(1 - 2b/3)} 10^{9.1 + 3M_g/2}, \quad (5)$$

where n is the number of terrestrial years, M_g is the magnitude of the seismic event and Γ is the gamma function. The geodetic moment rate, which depends on the strain rate originating from lithospheric deformation, ˙e, scales as

$$\dot{M}_{\text{geodetic}} = 2\dot{\varepsilon}\alpha\mu V. \quad (6)$$

In this equation, α is a seismic efficiency factor ranging from 0 to 1, V is the seismogenic volume and μ is the average shear modulus throughout this volume.

Given the non-observation of 3.9 magnitude marsquakes at high latitudes over the 4-year lifetime of the InSight mission, substituting Ṁ_geodetic and Ṁ_seismic in the above equations can be used to infer a maximum strain rate as a function of α, μ and V. The shear modulus μ is obtained from the interior models and the seismogenic volume V is calculated assuming a seismogenic depth given by the 573 or 1,073 K mantle isotherms⁵². In our models, μ and V are found to range from 44–72 GPa and 67–334 km, with a low shear modulus corresponding to models with a thin seismogenic layer thickness. The seismic efficiency factor, α, however, is mostly unknown for Mars, although it has been shown to be less than 0.7 (ref. 10). As a result, we consider this parameter to be 0.7 or 0.03, the latter of which is the minimum value that has been estimated on Earth⁴¹.

We derive the maximum strain rate from the non-observation of marsquakes in the northern regions over 4 years of InSight data collection to be 1.84 × 10⁻¹⁸ s⁻¹. As shown in Fig. 3, young polar caps (less than 1 Ma) show strain rates that are higher and inconsistent with the InSight-derived moment rate, and these models are thus excluded.

## Data availability

Interior models are available at Zenodo (https://doi.org/10.5281/zenodo.13710168)⁵³. Radar data are available from ref. 6.

## Code availability

ALMA is available at https://github.com/danielemelini/ALMA3 and Perple_X is available at https://www.perplex.ethz.ch/. The GAIA code is a proprietary code of the German Aerospace Center (DLR). Users interested in working with it should contact A.-C.P. (ana.plesa@dlr.de) and C. Hüttig (christian.huettig@dlr.de).

45. Boynton, W. V. et al. Concentration of H, Si, Cl, K, Fe, and Th in the low- and mid-latitude regions of Mars. *J. Geophys. Res. Planets* https://doi.org/10.1029/2007JE002887 (2007).

46. Ruedas, T. & Breuer, D. On the relative importance of thermal and chemical buoyancy in regular and impact-induced melting in a Mars-like planet. *J. Geophys. Res. Planets* **122**, 1554–1579 (2017).

47. McNutt, M. K. Lithospheric flexure and thermal anomalies. *J. Geophys. Res. Solid Earth* **89**, 180–194 (1984).

48. Broquet, A. Te_HF_Conversion: (0.2.3). Zenodo https://doi.org/10.5281/zenodo.6050262 (2022).

49. Connolly, J. A. D. The geodynamic equation of state: what and how. *Geochem. Geophys. Geosyst*. https://doi.org/10.1029/2009GC002540 (2009).

50. Fanara, L., Gwinner, K., Hauber, E. & Oberst, J. Present-day erosion rate of north polar scarps on Mars due to active mass wasting. *Icarus* **342**, 113434 (2020).

51. Thomas, P. C. et al. Mass balance of Mars' residual south polar cap from CTX images and other data. *Icarus* **268**, 118–130 (2016).

52. Plesa, A.-C. et al. Present-day Mars' seismicity predicted from 3-D thermal evolution models of interior dynamics. *Geophys. Res. Lett.* **45**, 2580–2589 (2018).

53. Brouquet, A. Data for glacial isostatic adjustment reveals Mars' interior viscosity structure. *Zenodo* https://doi.org/10.5281/zenodo.13710168 (2024).

54. Rudolph, M. L., Lekic, V. & Lithgow-Bertelloni, C. Viscosity jump in Earth's mid-mantle. *Science* **350**, 1349–1352 (2015).

**Acknowledgements** A.B. is funded by the Alexander von Humboldt Foundation. A.B., A.-C.P. and V.K. thank Geo.X—The Research Network for Geosciences in Berlin and Potsdam for providing financial support for the ‘Grow Your Idea!’ workshops.
**Author contributions** A.B. conceptualized the work and methodology, carried out the data analyses and modelling. A.-C.P. ran the thermal evolution simulations. V.K. and B.C.R. contributed to the viscoelastic deformation modelling. A.G. contributed to the reanalysis of the gravity signal associated with ice accumulation and isostatic adjustment rates. M.K. contributed to the seismic and geodetic moment analyses. All authors contributed to interpreting the results and writing the manuscript.

**Funding** Open access funding provided by Deutsches Zentrum für Luft- und Raumfahrt e.V. (DLR).

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-024-08565-9.

**Correspondence and requests for materials** should be addressed to A. Broquet.

**Peer review information** *Nature* thanks Isaac Smith and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_8.png" /><figcaption>**Extended Data Fig. 1** | Elastic lithosphere thicknesses of the 27 selected interior models. The grey shade indicates the allowed elastic thicknesses from radar analyses (Br20, ref. 6) and the black dots provide the range of possible</figcaption></figure>

elastic thicknesses. Line colours and interior models are the same as in Fig. 2.
Interior model names are described in the Methods and Extended Data Table 1.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Extended Data Fig. 2** | Present-day average interior structure of the northern regions (poleward of 60°N). Rigidity (a), viscosity (b), temperature (c), and density (d) of the interior as a function of depth. Line colours are the same as in Fig. 2 and Extended Data Fig. 1.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 3** | Effect of the ice-accumulation history on polar deformations. Time evolution of the maximum deformation beneath the north polar cap (a, b) considering different ice-accumulation histories (linear,

sinusoidal, sawtooth; c, d) for the 47-3.0_1.25λ (left) and 21-2.65 (right) models.
The accumulation histories in c and d are similar.
<figure><img src="image_15.png" /><figcaption>**Extended Data Fig. 4** | Time variable zonal degree 2 and 3 coefficients of the gravitational potential. Observed values and periodic fit (a, b). Residuals from the fit and linear trend in the residuals (c, d).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Extended Data Fig. 5 | Constrained interior temperature and viscosity.** Interior temperature (a) and viscosity (b) ranges of the northern regions (poleward of 60°N) and planetary average for accepted models.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

**Extended Data Fig. 6** | Effect of a periodic ice-accumulation history on polar deformations. Time evolution of the maximum deformation beneath the north polar cap (a, b) considering a 10 Myr periodic polar cap loading (c, d) for

the 47-3.0 1.25λ (left) and 21-2.65 (right) models. The accumulation histories in c and d are similar. The blue line provides a reference model with ice cap age of 10 Ma that neglects the past ice history.
<figure><img src="image_21.png" /><figcaption>**Extended Data Fig. 7** | Effect of the ice-accumulation history on polar deformations. Time evolution of the maximum deformation beneath the north polar cap (a, b) considering the loading history of the standard model in ref. 23 (c, d) for the 47-3.0_1.25λ (left) and 21-2.65 (right) models. The accumulation</figcaption></figure>

histories in c and d are similar. In this model, the past ice cap that decayed at ~8 Ma is assumed to have been long-lasting and stable over the last 500 Ma. The dashed line provides a reference model that neglects this past ice cap.
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

**Extended Data Fig. 8** | Effect of the ice-accumulation history on polar deformations. Time evolution of the maximum deformation beneath the north polar cap (a, b) considering the loading history of the standard model in ref. 23 (c, d) for the 47-3.0_1.25λ (left) and 21-2.65 (right) models. The accumulation

histories in c and d are similar. In this model, the past ice cap that decayed at ~8 Ma is assumed to have formed at 14 Ma. The dashed line provides a reference model that neglects this past ice cap.
<figure><img src="image_26.png" /><figcaption>**Extended Data Fig. 9 | Polar cap thickness and basement shape.** Thickness of the polar cap assuming a flat basement (a) and deformed basement of the polar cap assuming a real dielectric constant of 2.75 in our radar measurements (b).</figcaption></figure>
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

**Extended Data Fig. 10** | Love numbers for Maxwell and Andrade rheologies.
Mars interior response for Maxwell and Andrade rheology for the degree 8 (a) and 16 (b) load Love number $h'$ and associated interior viscosity (c) for two

tested models. Love numbers in this figure are calculated assuming a Heaviside load history. The Andrade rheology is defined following ref. 8 and assumes an Andrade creep parameter of 0.3.
<table><caption>Extended Data Table 1 | Summary of tested models. If not indicated, nominal parameters are initial temperature of $T_{init}=1650$ K, activation energy of $E=300$ kJ mol$^{-1}$, temperature difference across the mantle of $\Delta T=2000$ K</caption><thead><tr><th>Interior Model</th><th>ρ<sub>south</sub>, g cm<sup>-3</sup></th><th>ρ<sub>north</sub>, g cm<sup>-3</sup></th><th>d<sub>c</sub>, km</th><th>d<sub>c</sub><sup>InSight</sup>, km</th><th>λ</th><th>η<sup>500</sup>, Pa s</th><th>dT dz<sup>-1</sup>, K km<sup>-1</sup></th></tr></thead><tbody><tr><td>21-2.55</td><td>2.55</td><td>2.55</td><td>30</td><td>21</td><td>1.00</td><td>2.0e+21</td><td>3.2</td></tr><tr><td>21-2.55_1.75λ</td><td>2.55</td><td>2.55</td><td>30</td><td>21</td><td>1.75</td><td>3.3e+21</td><td>3.1</td></tr><tr><td>21-2.65</td><td>2.65</td><td>2.65</td><td>32</td><td>21</td><td>1.00</td><td>2.3e+21</td><td>3.2</td></tr><tr><td>21-2.65_1.75λ</td><td>2.65</td><td>2.65</td><td>32</td><td>21</td><td>1.75</td><td>4.1e+21</td><td>3.1</td></tr><tr><td>21-2.55-2.65</td><td>2.55</td><td>2.65</td><td>29</td><td>21</td><td>1.00</td><td>2.4e+21</td><td>3.1</td></tr><tr><td>21-2.55-2.65_1.75λ</td><td>2.55</td><td>2.65</td><td>29</td><td>21</td><td>1.75</td><td>2.9e+21</td><td>3.2</td></tr><tr><td>21-2.55-2.65_2λ</td><td>2.55</td><td>2.65</td><td>29</td><td>21</td><td>2.00</td><td>3.8e+21</td><td>3.1</td></tr><tr><td>31-2.55_2.1λ</td><td>2.55</td><td>2.55</td><td>41</td><td>31</td><td>2.10</td><td>1.5e+22</td><td>3.0</td></tr><tr><td>41-3.1_T<sub>init</sub></td><td>3.10</td><td>3.10</td><td>62</td><td>41</td><td>1.00</td><td>1.3e+22</td><td>2.5</td></tr><tr><td>41-3.1_E</td><td>3.10</td><td>3.10</td><td>62</td><td>41</td><td>1.00</td><td>1.2e+22</td><td>2.4</td></tr><tr><td>41-3.1_T<sub>init</sub>_E</td><td>3.10</td><td>3.10</td><td>62</td><td>41</td><td>1.00</td><td>8.3e+21</td><td>2.5</td></tr><tr><td>41-3.1_η<sub>jump</sub></td><td>3.10</td><td>3.10</td><td>62</td><td>41</td><td>1.00</td><td>3.0e+22</td><td>2.5</td></tr><tr><td>41-3.1_T<sub>init</sub>_η<sub>jump</sub></td><td>3.10</td><td>3.10</td><td>62</td><td>41</td><td>1.00</td><td>2.2e+22</td><td>2.6</td></tr><tr><td>41-3.1_E_ΔT</td><td>3.10</td><td>3.10</td><td>62</td><td>41</td><td>1.00</td><td>6.0e+21</td><td>2.6</td></tr><tr><td>47-2.6</td><td>2.60</td><td>2.60</td><td>58</td><td>47</td><td>1.00</td><td>4.0e+21</td><td>2.9</td></tr><tr><td>47-2.6-2.8</td><td>2.60</td><td>2.80</td><td>50</td><td>47</td><td>1.00</td><td>3.6e+21</td><td>3.0</td></tr><tr><td>47-2.8</td><td>2.80</td><td>2.80</td><td>61</td><td>47</td><td>1.25</td><td>1.9e+22</td><td>2.9</td></tr><tr><td>47-2.8-3.0</td><td>2.80</td><td>3.00</td><td>52</td><td>47</td><td>1.00</td><td>4.1e+21</td><td>3.0</td></tr><tr><td>47-3.0</td><td>3.00</td><td>3.00</td><td>69</td><td>47</td><td>1.00</td><td>8.0e+21</td><td>2.7</td></tr><tr><td>47-2.6-3.0</td><td>2.60</td><td>3.00</td><td>43</td><td>47</td><td>1.00</td><td>3.3e+21</td><td>3.0</td></tr><tr><td>47-2.6-3.0_1.38λ</td><td>2.60</td><td>3.00</td><td>43</td><td>47</td><td>1.38</td><td>1.8e+22</td><td>2.9</td></tr><tr><td>47-2.6-3.0_1.75λ</td><td>2.60</td><td>3.00</td><td>43</td><td>47</td><td>1.75</td><td>1.0e+22</td><td>2.6</td></tr><tr><td>47-2.6-3.0_2λ</td><td>2.60</td><td>3.00</td><td>43</td><td>47</td><td>2.00</td><td>2.5e+22</td><td>2.4</td></tr><tr><td>47-2.8-3.0_1.5λ</td><td>2.80</td><td>3.00</td><td>52</td><td>47</td><td>1.50</td><td>1.2e+22</td><td>2.5</td></tr><tr><td>47-2.8-3.0_1.75λ</td><td>2.80</td><td>3.00</td><td>52</td><td>47</td><td>1.75</td><td>6.4e+22</td><td>2.3</td></tr><tr><td>47-3.0_1.25λ</td><td>3.00</td><td>3.00</td><td>69</td><td>47</td><td>1.25</td><td>4.9e+22</td><td>2.3</td></tr><tr><td>49-3.0_0.83λ</td><td>3.00</td><td>3.00</td><td>71</td><td>49</td><td>0.83</td><td>1.8e+22</td><td>2.9</td></tr></tbody></table>

If noted in the interior model's name, $E=325$ kJ mol$^{-1}$, $T_{init}=1850$ K, and $\Delta T=2205$ K. The two models annotated with $\eta_{jump}$ have a 25-fold viscosity jump at around 900 km as suggested on Earth$^{54}$. The heat production factor ($\lambda$) scales the crustal heat production with respect to the nominal average Gamma ray estimate of 49 pW kg$^{-1}$ (ref. 28). $d_c$ is the average crustal thickness, $d_c^{InSight}$ is the crustal thickness at the InSight landing site, $\eta^{500}$ is the volume average present-day viscosity from 500 km to the core, and $dT dz^{-1}$ gives the upper mantle (100–400 km) thermal gradient. The three models that fit our constraints on glacial isostatic adjustment are displayed in bold.
Extended Data Table 2 | Summary of time-variable gravity analyses for the C<sub>20</sub> and C<sub>30</sub> coefficients

<table><thead><tr><th>Coefficient</th><th>m (×10<sup>-18</sup>)</th><th>b (×10<sup>-10</sup>)</th><th>p-value</th><th>r-value</th></tr></thead><tbody><tr><td>C<sub>20</sub></td><td>1.5±1.6</td><td>6.2</td><td>0.35</td><td>0.13</td></tr><tr><td>C<sub>30</sub></td><td>0.2±0.3</td><td>1.6</td><td>0.46</td><td>0.06</td></tr></tbody></table>

Values for *m* and *b* are for the *mx+b* linear fit to the residuals. The *p*-value tests the null hypothesis that the residuals have no correlation with time based on Wald's test, and the *r*-value provides the Pearson correlation coefficient of the residuals with time.