# Thickness of the Seasonal Deposits at the Martian North Polar Region From Shadow Variations of Fallen Ice Blocks

Haifeng Xiao, Yuchi Xiao, Shu Su, Frédéric Schmidt, Luisa M. Lara, Pedro J. Gutierrez

► To cite this version:

Haifeng Xiao, Yuchi Xiao, Shu Su, Frédéric Schmidt, Luisa M. Lara, et al.. Thickness of the Seasonal Deposits at the Martian North Polar Region From Shadow Variations of Fallen Ice Blocks. Journal of Geophysical Research. Planets, 2024, 129, 6, pp. 247-263. 10.1029/2023JE008078. insu-04843763

HAL Id: insu-04843763

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

Submitted on 17 Dec 2024

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
# JGR Planets

## RESEARCH ARTICLE

10.1029/2023JE008078

**Key Points:**

* We propose to examine the shadow variations of the ice blocks at the Martian polar region to infer the thickness of the seasonal deposits

* Maximum thickness of the seasonal deposits at the study scarp in MY31 is 1.63 ± 0.22 m to which snowfalls contribute 0.97 ± 0.13 m

* Seasonal deposits at the study scarp are up to 0.8 m shallower than previous measurements during spring

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

H. Xiao and Y. Xiao,

hxiao@iaa.es;

yuchi_xiao@cug.edu.cn

**Citation:**

Xiao, H., Xiao, Y., Su, S., Schmidt, F., Lara, L. M., & Gutierrez, P. J. (2024). Thickness of the seasonal deposits at the Martian north polar region from shadow variations of fallen ice blocks. *Journal of Geophysical Research: Planets*, 129, e2023JE008078. https://doi.org/10.1029/2023JE008078

Received 1 SEP 2023

Accepted 20 MAR 2024

**Author Contributions:**

**Conceptualization:** Haifeng Xiao

**Data curation:** Haifeng Xiao

**Formal analysis:** Haifeng Xiao, Frédéric Schmidt, Luisa M. Lara, Pedro J. Gutierrez

**Funding acquisition:** Haifeng Xiao, Luisa M. Lara, Pedro J. Gutierrez

**Investigation:** Haifeng Xiao, Yuchi Xiao

**Methodology:** Haifeng Xiao, Shu Su, Frédéric Schmidt, Pedro J. Gutierrez

**Project administration:** Haifeng Xiao, Luisa M. Lara, Pedro J. Gutierrez

**Resources:** Haifeng Xiao

**Software:** Haifeng Xiao, Shu Su

© 2024. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

# Thickness of the Seasonal Deposits at the Martian North Polar Region From Shadow Variations of Fallen Ice Blocks

Haifeng Xiao¹,² ID, Yuchi Xiao³ ID, Shu Su², Frédéric Schmidt⁴,⁵ ID, Luisa M. Lara¹, and Pedro J. Gutierrez¹

¹Instituto de Astrofísica de Andalucía (IAA-CSIC), Granada, Spain, ²Institute of Geodesy and Geoinformation Science, Technische Universität Berlin, Berlin, Germany, ³School of Mechanical Engineering and Electronic Information, China University of Geosciences (CUG), Wuhan, China, ⁴Université Paris-Saclay, CNRS, GEOPS, Orsay, France, ⁵Institut Universitaire de France (IUF), Paris, France

**Abstract** The seasonal deposition and sublimation of CO₂ constitute a major element in the Martian volatile cycle. Here, we propose to use the shadow variations of the ice blocks at the foot of the steep scarps of the North Polar Layered Deposits (NPLD) to infer the vertical evolution of the seasonal deposits. We conduct an experiment at a steep scarp centered at (85.0°N, 151.5°E). We assume that no snowfall remains on top of the selected ice blocks, the frost ice layer is homogeneous around the ice blocks and their surroundings, and no significant moating is present. We show that the average thickness of the seasonal deposits due to snowfalls in Mars Year 31 is 0.97 ± 0.13 m at Ls = 350.7° in late winter. The large depth measured makes us wonder if snowfalls are more frequent and violent than previously thought. Meanwhile, we show that the average frost thickness in Mars Year 31 reaches 0.64 ± 0.18 m at Ls = 350.7° in late winter. Combined, the total thickness of the seasonal cover in Mars Year 31 reaches 1.63 ± 0.22 m at Ls = 350.7° in late winter, continuously decreases to 0.45 ± 0.06 m at Ls = 42.8° in middle spring and 0.06 ± 0.05 m at Ls = 69.6° in late spring. These estimates are up to 0.8 m lower than the existing Mars Orbiter Laser Altimeter results during the spring. Meanwhile, we observe that snow in the very early spring of Mars Year 36 can be 0.36 ± 0.13 m thicker than that in Mars Year 31. This study demonstrates the dynamics of the Martian climate and emphasizes the importance of its long-term monitoring.

**Plain Language Summary** Like Earth, Mars also has seasons. Up to one third of the atmospheric CO₂ annually exchanges with the polar surface through seasonal deposition/sublimation processes. Deposition can be either atmospheric precipitation as snowfall or direct surface condensation as frost. At the steep scarps of the North Polar Layered Deposits (NPLD), fractured ice fragments can detach and fall to form ice blocks. We propose to use variations in the shadows of these ice blocks, observed in the High Resolution Imaging Science Experiment images, to infer the thickness evolution of the seasonal deposits. We make reasonable assumptions about the distribution of snowfall and frost around the ice blocks and their surroundings, which allow us to separately measure the thickness of snowfall and frost. Meanwhile, we introduce a novel approach that allows us to estimate the thickness of the seasonal deposits during late winter and early spring when image quality is insufficient. This approach also enables us to peer into the interannual thickness variations of snowfall. We carry out a successful experiment at a scarp centered at (85.0°N, 151.5°E). The obtained thickness measurements demonstrate the dynamics of the Martian volatile cycling and can be used to constrain the Martian climate models.

## 1. Introduction

When life emerged on Earth billions of years ago, the climate of Mars could have been with a thick atmosphere and a circumpolar ocean of liquid water in the northern hemisphere (Head III et al., 1999; Schmidt et al., 2022). Unfortunately, the present Martian surface is barren and arid with water tied up in the polar caps and in the subsurface, or even in the Martian crust (Scheller et al., 2021). The Martian Polar Layered Deposits (PLD) are predominantly pure water ice accumulated due to periodic orbital forcing (Becerra et al., 2017; Laskar et al., 2002), with a diameter of 1,000 km across and a thickness of a few kilometers (Nerozzi & Holt, 2019; Plaut et al., 2007). The spiral troughs that dissect the deposits contain thousands of visible ice layers with varying ratio of dust that record the Martian climate history in the Late Amazonian (Becerra et al., 2017; Grima et al., 2009). The carving of the troughs is probably related to prevailing katabatic winds that spiral due to the Coriolis effects
**Supervision:** Haifeng Xiao,
Frédéric Schmidt, Luisa M. Lara, Pedro
J. Gutierrez
**Validation:** Haifeng Xiao
**Visualization:** Haifeng Xiao, Yuchi Xiao
**Writing – original draft:** Haifeng Xiao
**Writing – review & editing:**
Haifeng Xiao, Shu Su, Frédéric Schmidt,
Luisa M. Lara, Pedro J. Gutierrez

(I. B. Smith & Holt, 2010). The PLD, together with the outlying crater ice deposits (McGlasson et al., 2023; Sori et al., 2022), are the most accessible planetary climate records outside Earth, making them among the most compelling science targets in the Solar System (Becerra et al., 2021; I. B. Smith et al., 2020). On the steep scarps at the margins of the North Polar Layered Deposits (NPLD), mass-wasting processes like dust-ice avalanches and ice block falls have been frequently observed (Fanara et al., 2020b; Russell et al., 2008; Su, Fanara, Xiao, et al., 2023; Su, Fanara, Zhang, et al., 2023). The latter is in the form of scattered ice blocks and accumulating apron-like debris at the foot of the scarps. These phenomena represent multiple alternative modes of erosion in addition to sublimation, which are held accountable for mass losses of the polar ice caps. Mass loss rate due to these phenomena can compete with the outward motion due to viscous deformation to shape the NPLD's rheology and evolution (Fanara et al., 2020b; Sori et al., 2015). In terms of the ice blocks, ice fragments that fall from the steep scarps are mainly protrusions that break due to structural failure once being subject to internal/external triggers, for example, thermal expansion/contraction and then stress-induced polygonal fracturing, loading/ removal of the topping seasonal deposits (Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022), or even downward katabatic winds (Byrne et al., 2017).

Due to its obliquity (currently at ~25° compared to 23.5° of Earth), Mars experiences seasons as the planet orbits the Sun. When the temperature drops below the CO₂ condensation point (approximately 148 K at average Martian surface pressure, but can range from 130 to 154 K depending on elevation and season) in the fall and winter, the CO₂ solidifies and accumulates as snow or frosts in the polar regions. Then when the temperature increases during spring, these deposits sublimate into the atmosphere. Each Martian year, up to one third of the CO₂ in the atmosphere gets involved in the seasonal deposition/sublimation process (Leighton & Murray, 1966). Two depositional mechanisms exist, that is, atmospheric precipitation as snowfall and direct surface condensation as frost (Määttänen & Montmessin, 2021). The resultant Seasonal South/North Polar Caps (SSPC, SNPC) at their maximum can reach around 50°S/N in the beginning of winter under current conditions (Piqueux et al., 2015). Understanding these seasonal processes can place vital constraints on the Martian volatile cycles, and help with the determination of the current mass balance of the polar ice reservoirs (Becerra et al., 2021). The thickness of the surficial layer would have important implications for the feasibility or trafficability of future landers, rovers, or helicopters that would drill into the PLD and decipher the stored paleoclimate of Mars (Matthies et al., 2022; I. B. Smith et al., 2020). Through Mars Climate Sounder (MCS) observations, it is now believed that water ice particles in the atmosphere can act as condensation nuclei onto which the CO₂ ice condenses and they can be deposited onto the polar caps along with the CO₂ snowfalls (Alsaeed & Hayne, 2022). Thus, by constraining the quantities of the CO₂ snowfalls, we can also gain insights into the amount of water that can be annually removed from the atmosphere through this scavenging process. This information can also help to determine the current mass balance of the residual north polar cap, that is, whether it is accumulating or ablating. Dynamic geological phenomena formed during sublimation of the SSPC/SNPC, for example, dark fans, polygonal cracks, spiders (south polar region)/furrows (north polar region), and alcoves in the dune fields, can be better modeled and interpreted given meaningful thickness evolution measurements of the overlying seasonal layer (e.g., Dundas et al., 2021; Hansen et al., 2013; Mc Keown et al., 2023; Portyankina et al., 2010; Schmidt & Portyankina, 2018). Meanwhile, the annual growth and retreat of the seasonal deposits are capable of inducing significant elastic displacements of the lithosphere (Métivier et al., 2008). Quantifying these displacements through future dedicated geodetic missions can place important constraints on the current thermal and rheologic state of the Martian interior (Wagner et al., 2023). The changing mass loads can also cause small but measurable effects on Mars gravity and—through changing mass distribution and moments of inertia—rotation (i.e., nutation, polar motion, and length of day variation; Defraigne et al., 2000; Le Maistre et al., 2023; Van den Acker et al., 2002).

The direct depth variation measurements of the seasonal deposits have been made by exploiting dynamic Mars Orbiter Laser Altimeter (MOLA) elevation profiles (Aharonson et al., 2004; D. E. Smith et al., 2001; Xiao, Stark, Schmidt, Hao, Su, et al., 2022; Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022). The MOLA estimates can be easily extended to cover the entire polar regions. D. E. Smith et al. (2001) calculated the height differences of the MOLA heights to median-filtered reference surface at various latitudinal annuli and measured the maximum thickness to be ~1 m at both poles. Aharonson et al. (2004) fitted sinusoidal curves to the height differences obtained at locations where two MOLA profiles intersect, that is, cross-overs, and estimated the maximum depth variations to be ~1.5 m at the north polar region and ~2.5 m at the south polar region. Recently, by reprocessing the MOLA profiles and self-registering them (Xiao, Stark, Steinbrügge, et al., 2021; Xiao, Stark, Steinbrügge, et al., 2022). Xiao, Stark, Schmidt, Hao, Su, et al. (2022) and Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2022).
derived both spatial and temporal thickness variations of the seasonal polar caps with a maximum of about 2.5 m at the south and 1.3 m at the north. In particular, they brought attention to abnormal thickness patterns of the SNPC at Olympia Undae, where peak-to-peak thickness variations up to 4 m and off-season thickness increases and decreases up to 3 m in magnitude have been observed. Olympia Undae, which roughly spreads from 78°N to 83°N in latitude and 120°E to 240°E in longitude (refer to Figure 4 for its location), is the largest continuous dune field on Mars (Hayward et al., 2010). It is a part of the polar erg that denotes the vast circumpolar dark dune fields. However, it should be noted that the MOLA altimetric dataset is limited in time to MY24 and MY25, preventing it to resolve inter-annual changes in the seasonal snow/ice depth.

Andrieu et al. (2018) applied Bayesian inversion techniques involving a radiative transfer model to constrain the impurity content and depth of the SSPC at a dune field of Richardson Crater (72°S, 180°E), and a maximum thickness of ~0.4 m was found in early southern spring. Currently, they are trying to extend the measurements to the SSPC over the “cryptic region” that features translucent CO₂ slab ice and cold-jetting during the southern spring (refer to Hansen et al. (2010) for the location of this enigmatic region). Raguso and Nunes (2021) performed advanced radar processing techniques to the SHAllow RADar (SHARAD) dataset to garner the best possible resolution and Signal-to-Noise Ratio (SNR) and co-registered the subsurface reflectors for estimation of the two-way time delay differences, hence the thickness of the Martian seasonal layer. Unlike MOLA, these measurements were not affected by ephemeris errors of the spacecraft but can be biased by possible presence of slope differential between the surface and the subsurface reference reflectors. Unfortunately, the experimental results are currently unavailable and extension to other regions with rougher surfaces and less distinct subsurface reflectors can be difficult.

Another viable way is to relate the rock shadow length changes in high-resolution optical images to that of the snow/ice depth. The High Resolution Imaging Science Experiment (HiRISE) camera on-board the NASA's Mars Reconnaissance Orbiter (MRO) operates in visible wavelengths, and with a telescopic lens that produces images with a pixel size of 0.25–1.3 m (McEwen et al., 2007). This high-resolution can enable accurate identification of rocks as small as 0.7 m in diameter. The spacecraft operates in a near sun-synchronous orbit, providing images of similar solar azimuth, beneficial for comparison of the shadow lengths that indicates snow/ice depth variations. A HiRISE Digital Elevation Model (DEM) can be derived from a geometric stereo pair acquired on different orbits so that a moderate convergence angle between the two viewing directions is formed (~10° to ~25°, McEwen et al., 2007). The camera started to continuously acquire high-resolution pictures of targeted regions on 29 September 2006, obtaining 5–20 observations per day, which has led to a high-cadence time series of imagery. Cull et al. (2010) utilized HiRISE images to show that the CO₂ snow/ice thickness reduced from ~30 cm in early spring to less than 5 cm by middle spring at the Mars Phoenix landing site (68.22°N, 125.70°W). Mount and Titus (2015) also used rock shadow measurements in HiRISE images to infer the seasonal snow/ice depth at three sites of distinct morphologies (at latitudes between 68°N and 75°N). They showed that the effects of moats (sublimation of snow/ice proximal to rocks) and crowns (accumulation of snow/ice on the tops of rocks) can significantly modify rock shadow measurements and hence snow/ice depth (refer to Figure 3 of Mount and Titus (2015) for a schematic). Therefore, they applied an empirical model for correction and performed error propagation analysis to reflect the uncertainty in the measurements. These snow/ice depth measurements were then combined with visible and thermal observations to calculate the bulk density of the seasonal ice cover over time. Despite its high precision, this approach can be spatially limited to fields where rocks are present on the surface. Unfortunately, at high polar latitudes where maximum snow/ice accumulation happens, rocks are unavailable for the purposes of looking into the seasonal snow/ice thickness evolution.

As aforementioned, the only existing measurement of the seasonal snow/ice depth variations at high latitudes come from MOLA records, which date back to MY24/25. Here, we propose to use shadows of the ice blocks at the foot of the steep scarps as observed in the HiRISE images, complementing that of the rocks, as an alternative way to infer the depth of the seasonal deposits. We show how to relate the length of the ice block shadows to its height using a rigorous geometric model, which is based on orthorectified images and takes both the solar and surface properties into consideration. Building on this model, we present two independent and complementary approaches to shed light on the thickness evolution of the ephemeral deposits: (a) Thickness from Height Variations (THV) that subtracts the ice block heights measured in the summer when free of seasonal cover to that in the spring; (b) Thickness from Height Bounds (THB) that locates ice blocks that have been completely covered to place lower bounds on the thickness of the seasonal deposits, and ice blocks that have not been completely submerged to put upper bounds. We carry out the experiments at a scarp centered at (85.0°N, 151.5°E) and show
the feasibility of these applications. We note that while THV is temporally limited to mid-to-late spring, THB is
capable of yielding measurements in late winter and early spring. Starting from late northern summer, thin haze
can rapidly grow into thick water ice clouds, the so-called polar hood. This large-scale feature can persist from
late summer, fall, and all the way to winter, and even early spring (Benson et al., 2011; Brown et al., 2016; Calvin
et al., 2015; Navarro et al., 2014). Fortunately, plenty of unobscured HiRISE images are available during late
winter and early spring. During experiments with these proposed methods, we observe that moats do not exist
around the ice blocks and the depth of the crowns over the ice blocks quasi-linearly decreases to zero when the
seasonal deposits completely sublimate back into the atmosphere. Thus, the empirical correction scheme
described in Mount and Titus (2015), that by adding the snow/ice thickness increases between consecutive
seasons to all prior-season thicknesses, can be inapplicable in our case. We hence make reasonable assumptions
and propose to use the widening of the ice blocks as a proxy to approximate and correct for the depth of the crowns
over the ice blocks. These assumptions also enable us to decompose the contributions of the snowfalls and direct
condensation to the thickness of the seasonal ice deposits and estimate them separately. Our ultimate goal is to
apply these approaches to all active scarps at high polar latitudes, and rock fields at lower polar latitudes, to obtain
good samplings of the vertical evolution of the SNPC. The expected results can also serve to calibrate existing
MOLA results and validate contemporary anticipated SHARAD results.

The paper is structured as follows: In Section 2, bundle adjustment and orthorectification of the HiRISE images are introduced. Based on the orthorectified images, the shadowing model of ice blocks resting at the foot of the steep scarps is presented, and the THV and THB approaches to obtain thickness evolution of the SNPC are described. The study scarp centered at (85.0°N, 151.5°E) is introduced in Section 3. These are followed by the application of the two independent approaches to ice blocks at the study scarp (Section 4). After that, the precision of the THV results, the automation of the THB approach, the interannual variations of snowfalls, the estimation and correction for direct condensation effects, the possible biases in the MOLA-derived thicknesses, the comparison of measured snowfall thickness to that predicted by a snowing model, and the prospects of future work are successively presented in Section 5. Finally, conclusions are drawn in Section 6.

**2. Methods**

To correctly relate an ice block's shadow length to its height, bundle adjustment and subsequent orthorectification of the images should be implemented to remove the image distortions due to oblique viewing angle and topographic relief. For the orthorectification to be carried out, a precise DEM has to be made using image matching of the bundle-adjusted stereo pairs. Previous studies using the objects' shadows to infer their heights focused on relatively flat regions without significant undulations and thus assumed the surface to be a horizontal plane (e.g., Blackburn et al., 2010; Cull et al., 2010; Mount & Titus, 2015; P. C. Thomas et al., 2016). However, in our case, the slopes in the regions where ice blocks reside in the Basal Unit outcrops, immediately underlying the NPLD, are significant (up to ~30°) and have to be considered in the shadowing geometry. The established shadowing model is then used in the proposed THV and THB approaches to determine the thickness evolution of the seasonal deposits.

**2.1. Bundle Adjustment, DEM Generation, and Image Orthorectification**

We use the raw and unprocessed HiRISE images, that is, the Experimental Data Record products, to do bundle
adjustment and produce a DEM of the study region using a stereo pair and the Ames Stereo Pipeline software
(Beyer et al., 2018; Hepburn et al., 2019). This DEM is then used to orthorectify all the available images. We note
that Reduced Data Record (RDR) products also exist that are radiometrically corrected images resampled to a
standard map projection. However, these RDR images lack the required geometric stability for stereo processing.
The sequential procedures applied are: (a) Mosaicking the individual 10 Charged-Coupled Devices (CCDs)
together to single images which includes de-jittering and radio-calibration; (b) Bundle adjusting all available
images to correct for errors in camera position and orientation (extrinsics only) and make them internally
consistent. Feature points are matched across images. A feature point can be identified in multiple overlapping
images, it is equivalent to the concept that a bundle of light rays must intersect at a single triangulated point on the
ground. In reality, the intersection can be imperfect due to residual errors. In bundle adjustment, each triangulated
ground point is projected back into the cameras. Then, the sum of squares of residuals between the pixel co-
ordinates of the feature points and the locations of the projected points, that is, reprojection errors, are minimized
through a robust least squares solver:
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Geometry schematic that illustrates how the solar azimuth, elevation, magnitude and aspect of the slope affect the shadow of the Ice Block situated perpendicular to the slope plane (O-ABC). Ph is the projected point of the Ice Block top onto the horizontal plane, while Ps is the top onto the inclined plane with a slope of β. The angle α denotes the elevation angle of the Sun. The angle ω represents the relative azimuth angle of the solar rays with respect to that of the slope.</figcaption></figure>

$$ \{\hat{\mathbf{R}}, \hat{\mathbf{t}}, \hat{\mathbf{P}}\} = \min \sum_{i=1}^{n} \sum_{j=1}^{m} (\mathbf{u}_{ij} - \pi(\mathbf{C}_j, \mathbf{R}_j \mathbf{P}_i + \mathbf{t}_j))^2, \quad (1) $$

where $\mathbf{u}_{ij}$ is the observed tie point image coordinate in pixel, $\mathbf{P}_i$ is the 3D coordinate of the ith ground point (a total of n), and $\mathbf{C}_j$ is the center coordinate of the jth camera (a total of m). Meanwhile, $\mathbf{R}_j$ and $\mathbf{t}_j$ denote the rotation and translation operations to the jth camera, $\pi(\mathbf{C}_j, \mathbf{R}_j \mathbf{P}_i + \mathbf{t}_j)$ is the reprojection operator that obtains the reprojected image coordinate in pixel of the ground point. Parameters to be adjusted are listed on the left of the equation; (3) Locating conjugate feature points through image matching techniques and derive corresponding disparity values. Sub-pixel correlation in image matching is performed to refine the disparity map, which is then converted to heights of the object ground points; (4) Gridding the point cloud of object heights to a DEM. This DEM is then applied to orthorectify all of the HiRISE images in the depositional area of the ice blocks. We note that the bundle adjustment lacks absolute ground control points and only tie points are used to improve the internal consistency of the images. As a result, there can exist lateral shifts between the DEMs generated in Su, Fanara, Zhang, et al. (2023) using stereo pairs from previous bundle-adjustment and the images from the bundle adjustment in this study. For more details of these processes, refer to the pipeline explained in Su, Fanara, Zhang, et al. (2023).

## 2.2. Relating the Height of an Ice Block to Its Shadow Length

Figure 1 shows the schematic on the shadowing geometry that we utilize to establish the relationship between the height of the Ice Block (H) and the measured length of the cast shadow as seen in an orthorectified HiRISE image (OPs'). The solar elevation angle is α, while the slope is β. We align the local coordinate system with the bearing of the sloped plane (O-ABC). The angle ω measures the angular separation between the sunlight and the orientation of the slope. Without considering the slope, the Ice Block in question casts a shadow OPh on the horizontal plane, coordinates of which can be written as
$$ \begin{cases} x_{Ph} = H \frac{\sin \omega}{\tan \alpha \cos \beta} \\ y_{Ph} = H \frac{\cos \omega}{\tan \alpha \cos \beta} \end{cases} \qquad (2) $$

Now, we consider the existence of the sloped plane and assume the intersection point of the light over the tip of the Ice Block with the sloped plane to be Ps. As Ps has to be somewhere between I and Ph, its coordinates can be like

$$ \begin{cases} x_{Ps} = t \times H \frac{\sin \omega}{\tan \alpha \cos \beta} \\ y_{Ps} = t \times H \frac{\cos \omega}{\tan \alpha \cos \beta} \\ z_{Ps} = (1-t) \times H \frac{1}{\cos \beta} \end{cases} \qquad (3) $$

where t is a scalar between 0 and 1. As can be seen in Figure 1, the ratio of z_Ps and y_Ps equals tan β, which leads to

$$ t = \frac{\tan \alpha}{\cos \omega \tan \beta + \tan \alpha} \qquad (4) $$

and the coordinates of the intersection point on the sloped plane can be written as

$$ \begin{cases} x_{Ps} = H \frac{\sin \omega}{\cos \omega \sin \beta + \tan \alpha \cos \beta} \\ y_{Ps} = H \frac{\cos \omega}{\cos \omega \sin \beta + \tan \alpha \cos \beta} \\ z_{Ps} = H \frac{\cos \omega \tan \beta}{\cos \omega \sin \beta + \tan \alpha \cos \beta} \end{cases} \qquad (5) $$

The length of the shadow on the slope as measured from the HiRISE orthorectified image is then

$$ L_S = \sqrt{x_{Ps}^2 + y_{Ps}^2} = \frac{H}{\cos \omega \sin \beta + \tan \alpha \cos \beta} \qquad (6) $$

Solar azimuth together with the elevation at the acquisition time of the HiRISE image can be calculated by exploiting the information stored in the Spacecraft, Planet, Instrument, Camera-matrix, and Events (SPICE) kernels (Acton, 1996). Specifically, we access pck00009.tpc Planetary Constant Kernel for orientation and shape of Mars. Meanwhile, we use the mar063.bsp file of the generic planet ephemeris to estimate relative position of Mars with respect to the Sun at a given time stamp. Slope magnitude and aspect at the Ice Block's location can be approximated by evaluating the elevation values within a 3 × 3 window in the HiRISE DEM generated in this study. The azimuthal difference of the solar illumination with respect to bearing of the local slope, ω that falls within [0°, 180°], can be related to the measured solar azimuth (φ) and slope aspect (μ) as shown in Figure S1 in Supporting Information S1. For the north polar stereographic projection centered at the North Pole, the projected body-fixed coordinates of the Ice Block with latitude φ and longitude L can be written as

$$ \begin{cases} x_{bf} = \frac{2R \cos \phi \sin L}{1 + \sin \phi} \\ y_{bf} = \frac{2R \cos \phi \cos L}{1 + \sin \phi} \end{cases} \qquad (7) $$

Once we know the map coordinates of the Ice Block, the analytical expression of ω can be derived as follows:
$$ 
\begin{aligned}
\omega &= |\pi/2 - (\varphi - \varepsilon - \pi/2) - \mu| = \left| \pi - \left(\varphi - \arctan \frac{x_{bf}}{y_{bf}}\right) - \mu \right| \\
&= |\pi - (\varphi - (\pi - L)) - \mu| \\
&= |2\pi - \varphi - L - \mu|,
\end{aligned} 
\quad (8) $$

where the angle $\varepsilon$ denotes the intersection angle between the map north and the direction for the Ice Block in question to the North Pole (Figure S1 in Supporting Information S1).

## 2.3. The THV Approach

By measuring the shadow lengths of the ice blocks in the orthorectified images, combined with auxiliary information on the solar and slope conditions, we can infer the ice block heights above the snow/ice cover during spring or the bare ground surface during summer using Equation 6. The HiRISE images can span eight Mars Years (MY), or more than 14 Earth years, during which the heights of the ice blocks can slowly shrink due to tumbling and ablation, for example, aeolian erosion and sublimation. Thus, we use the ice block heights obtained at the temporally adjacent northern summer as the reference to infer the depth of the seasonal cap in the springtime. When there exist multiple ice block height measurements within a single summer, their average is taken as the reference value, as a way to reduce measurement errors. See a detailed discussion in Section 5.1. Finally, acquired time-dependent depth values are plotted against the solar longitude for examination of their evolution patterns. Here, solar longitude (Ls) is used to express the seasonality on Mars and 0° < Ls < 90° stands for northern spring, 90° < Ls < 180° for northern summer, 180° < Ls < 270° for northern fall, and 270° < Ls < 360° for northern winter. We term this approach as Thickness from Height Variations (THV), as it basically subtracts the ice block heights measured in the summer to that measured in the spring to gauge the seasonal snow/ice depth in springtime.

Criteria for selection of the ice blocks as a reference to invert for the thickness variations of the seasonal deposits are as follows: (a) Ice blocks are high enough as not to be submerged by the seasonal snow/ice; (b) Ice blocks shaped like triangular prisms are the ideal candidates. The shadows in these cases run parallel to the ridge lines of the ice blocks (see Ice Block 3 in Figure 6 as an example). Spiky ice blocks, and so their shadows, are also considered solid options. The peaks on the ice blocks and in their shadows are normally easy to identify (examples are Ice Blocks 1 and 2 in Figure 6). In both of these cases, the effects of accumulation of snow particles over the ice block tops are minimized. However, it should be noted that the formation of a layer on the tops of the ice blocks, that is crowning, due to direct condensation of frosts can still take place (Section 5.4); (c) The ice blocks can tumble under gravity of its own and the seasonal layer, or even be triggered by the katabatic winds down from the polar caps. Furthermore, ablation like erosion due to winds can slowly shrink the size of the ice blocks. Thus, it is required that the morphology of the ice blocks and their heights should not significantly change during the bracketing summers. For triangular-prism-like ice blocks, the lengths of the shadow are treated as the distances between the parallel ice block crest lines and shadow edge lines along the solar azimuth at the acquisition times of the images (Figure S2 in Supporting Information S1). For peaked ice blocks, the lengths of the shadow are measured as that of the lines connecting tips of the shadows and their corresponding points along the terminators, that is, boundary lines separating the sunlit and shadowed portions of the ice blocks (Figure S2 in Supporting Information S1). When the latter is blurry, then we fix the starting points of the measuring lines at the shadow tips, fine-tune them to be parallel to the solar azimuth, and determine their ending points as the intersections between the lines and the respective terminators.

## 2.4. The THB Approach

Although HiRISE images without the impact of the polar hood can be acquired during late winter and early spring, the outlines of the shadows cannot be unambiguously identified due to unfavorable solar elevation angles, leading to failed application of the proposed THV approach. Meanwhile, the bundle adjustment can fail with images acquired in late winter or early spring as small-scale features stay covered by the thick seasonal deposits.

Here we devise an independent approach for measuring the depth of the seasonal deposits, which is capable of inferring the thickness during late winter and early spring (Figure 2). This approach takes advantage of relatively large quantity of ice blocks that cluster in specific regions (see the inset in Figure 5 for an example). The main idea
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Illustration of an alternative and independent approach employed to constrain the depth of the seasonal deposits (THB). Instances of ice blocks of various sizes and heights are randomly placed over the Basal Unit. It should be noted that this plot is a simplified illustration as in reality smaller ice blocks are much more frequently observed than larger ones (Fanara et al., 2020b). Meanwhile, only peak-shaped ice blocks are presented so as we can assume no snow particles on top of the ice blocks. $\Omega_f$ denotes the varying thickness of the seasonal layer at different epochs.</figcaption></figure>

behind it is that an ice block that has been completely covered means that the thickness of the seasonal layer exceeds the height of the ice block (a lower bound). Similarly, an ice block that has not been completely covered means that the thickness of the seasonal layer is lower than the ice block height (an upper bound). However, to place the strongest constraints, we have to locate the locally highest ice blocks that have been completely sub- merged and the locally lowest ice blocks that still stick out the frozen layer. We only select peaked and triangular- prism-shaped ice blocks to minimize the impact of accumulation of snow particles over the ice blocks which

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Schematic showing the offset inherited in the derived bounding constraints on the thickness of the seasonal layer at $L_s = \gamma^\circ$ in late winter or spring from the THB approach. IB 2 is an example of the locally smallest ice block visually determined to be uncovered by the seasonal layer to place an upper bound on its thickness ($\Omega_\gamma$). Meanwhile, IB 4 is an example of the locally largest ice block visually identified to be completely covered with the aim to place the corresponding lower bound. IBs 1 and 5 marked in red are representative ice blocks that are respectively too large and too small to put tight bounds on the thickness of the seasonal layer, which are thus not included in our analysis. The criterion by which we visually determine if an ice block is fully covered or not is if its shadow length is greater than the identification threshold of 2.83 pixels, or 0.71 m when the spatial resolution of the image is 0.25 m. Thus, the height of the bounding ice block during summertime ($H_m$, $m = 1, 2, 3, 4, 5$) always offsets the thickness of the seasonal deposits by a magnitude of $\xi_\gamma$. This offset can be shown in the ideal case of IB 3 where its shadow length exactly equals the threshold of 2.83 pixels, which can then precisely constrain the seasonal deposits' thickness to be $\Omega_\gamma = H_3 - \xi_\gamma$. The offset, that is $-\xi_\gamma$, can be related to the spatial resolution of the image and the solar elevation angle as in Equation 10. The boundaries for different types of ices are just for illustration purposes as there also exists directly condensed and CO$_2$-snowfall-scavenged water ice within the seasonal deposits (Alsaeed & Hayne, 2022; Appéré et al., 2011).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Map of the Martian north polar region with elevation represented by the reference Digital Elevation Model (DEM) from reprocessed and then self-registered Mars Orbiter Laser Altimeter (MOLA) altimetric profiles (Xiao, Stark, Schmidt, et al., 2021). The spatial resolution of the DEM is 1 km/pixel. The elevation is with respect to the Martian ellipsoid with an equatorial radius of 3,396.19 km and a mean polar radius of 3,376.20 km. The coverage gap poleward of 87°N is due to the absence of nadir-pointing profiles, a limitation related to MGS's orbital inclination. The map projection adopted is polar stereographic centered at the North Pole. The Residual North Polar Cap (green polygons), Chasma Boreale, and the study site (Scarp 1) are marked.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Illustration of the North Polar Layered Deposits (NPLD) and ice blocks at the foot of the steep scarp (Scarp 1) by a 3D view with the High Resolution Imaging Science Experiment (HiRISE) image acquired during the summertime of MY29 (PSP_009648_2650_RED) draped on the HiRISE Digital Elevation Model (DEM) generated in this study (no vertical exaggeration). Elevation difference from top of the NPLD to the Basal Unit at the bottom is ~600 m. Site is centered at (85.0° N, 151.5°E). Wide patches and aprons of debris are readily visible at the foot of the scarp. Inset map shows typical examples of ice blocks at foot of the scarp. It should be noted that the Basal Unit outcrops where ice blocks reside can still feature significant terrain slope (up to ~30°). We mark the locations of Ice Blocks 1, 2, and 3 used for the purpose of this study. For reference, the distance between Ice Blocks 1 and 3 is 4,564 m.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Images of Ice Blocks 1, 2, and 3 (centered at the sub-frames) during different seasons from late winter to summer. The images used are ESP_032856_2650_RED at Ls = 359.6° in MY31, ESP_024509_2650_RED at Ls = 17.3° in MY31, ESP_016439_2650_RED at Ls = 44.0° in MY30, ESP_016716_2650_RED at Ls = 53.6° in MY30, ESP_034610_2650_RED at Ls = 62.5° in MY32, and ESP_053730_2650_RED at Ls = 113.9° in MY34, respectively. The projection adopted is polar stereographic centered at the North Pole. The direction to the North Pole are marked. Illumination is to the bottom-left corner. The scale bar applies to all of the images. Note that ice blocks in springtime feature no circumferential shadows that are indicative of moating.</figcaption></figure>

otherwise can bias the thickness inversion (Section 5.4). As shown in Figure 2, snow/ice deposits are thickest in late winter and early spring with a thickness of $\Omega_1$, then as the insolation increases with time, $CO_2$ starts to sublimate and the thickness starts to decline and reach $\Omega_2$ at mid-spring and $\Omega_3$ in late spring, respectively. Surface is free of the seasonal layer during the adjacent summer where the heights of the ice blocks can be measured using the geometric model presented in Section 2.2. For $\Omega_1$, we can utilize Ice Blocks (IBs) 4 and 1 to place the lower and upper bounds, respectively. IB 4 is the locally highest ice block that has been completely covered during late winter and early spring, and IB 1 is the locally lowest ice block that has not been completely submerged. In a similar way, we can locate the ice blocks, marked using corresponding colors in Figure 2, to place bounds on $\Omega_2$ and $\Omega_3$, respectively. The mathematical expressions can be put as follows:

$$ \begin{cases} \Omega_1 \in [H_4, H_1] \\ \Omega_2 \in [H_5, H_6], \\ \Omega_3 \in [H_2, H_3] \end{cases} \qquad (9) $$

where $H_m$ denotes the height of IB $m$ measured in the adjacent summer. We term this alternative approach as Thickness from Height Bounds (THB) for simplicity and to distinguish it from the THV approach that subtracts the ice block heights in the summer to that in the spring (Sections 2.3).
Attention should be drawn to the fact that HiRISE limits the detection of shadows under $2\sqrt{2} = 2.83$ pixels in length (0.71 m for a pixel size of 0.25 m). This is a criterion tighter than the usually considered of 3 pixels as in reality the shadows can also misalign with the image lines and samples. For example, they can be right at the diagonals of the image grid. In practice, an ice block is visually determined to be fully covered as long as its shadow is less than 2.83 pixels in length. As such, when we determine a locally shortest ice block that has not been completely submerged or a locally tallest ice block that has been completely covered, it has already stuck out the seasonal layer at least or at most by

$$ \xi_Y = 2\sqrt{2} \times r_s \times \tan \alpha_Y. \qquad (10) $$

Here, $\xi_Y$ denotes magnitude of the offset at Ls = $\Upsilon^\circ$ during late winter or spring. The variable $r_s$ represents the spatial resolution of the images which is mostly 0.25 m but can also be 0.5 m for that acquired in late winter and early spring. The symbol $\alpha_Y$ denotes the solar elevation angle at the acquisition of the image with seasonal snow/ice cover. A schematic illustrating the offset is shown in Figure 3. We correct the derived thickness bounds for the offset of $-\xi_Y$ otherwise they would be systematically higher than that derived from the THV approach.

For illustration purposes, abstract probability distributions of bounding constraints at specific solar longitudes are computed using the Kernel Density Estimation (KDE; Silverman, 1986). In our case, we perform the KDE process using the Gaussian kernel ($\kappa$) and use the Silverman's rule to compute an optimal bandwidth of $\sigma$ to ensure moderate smoothness of the obtained density distribution:

$$ \rho_{\kappa}(x) = \frac{1}{N\sigma} \sum_{i=1}^{N} \kappa(x - x_i; \sigma) \quad \text{with} \quad \kappa(x; \sigma) \propto \exp\left(-\frac{x^2}{2\sigma^2}\right), \qquad (11) $$

where N is the total number of ice block height measurements, that is $x_i$, used to set upper or lower limits on the seasonal snow/ice depth. Violin plots are used for visualization purpose. We have also tried the Improved Sheather Jones algorithm for choosing the bandwidth which should be used when data is far from normal or multimodal (Botev et al., 2010). Unfortunately, the limited number of available constraints in our case prevents its effective application.

For statistically significant number of constraints, the thickness of the seasonal deposits at a specific solar longitude can be expected to fall within the interval formed by medians of the lower and upper bounds. To obtain more realistic uncertainty estimates, we take the standard errors of the median bounds into consideration using a scaled Median Absolute Deviation (SMAD; Leys et al., 2013). The SMAD is related to the median of the absolute deviations from the samples' median as

$$ \text{SMAD}_{\Upsilon} = f_s \times M[|\Omega_{\Upsilon} - M[\Omega_{\Upsilon}]|], \qquad (12) $$

where $\Omega_{\Upsilon}$ denotes the sample vector which contains all the thickness constraints at Ls = $\Upsilon^\circ$ in late winter or spring, $M[$] is the median operator, and $f_s = 1.4826$ is the scale factor. The scaled metric can be treated as a consistent estimator similar to the standard deviation of a Gaussian distribution. We extend the median lower bound downwards and the median upper bound upwards to obtain the adjusted bounding range:

$$ \Theta_{\Upsilon}^{\Omega} = \left[ M_{\Upsilon\_lb} - \frac{3 \times \text{SMAD}_{\Upsilon\_lb}}{\sqrt{N_{\Upsilon\_lb}}}, M_{\Upsilon\_ub} + \frac{3 \times \text{SMAD}_{\Upsilon\_ub}}{\sqrt{N_{\Upsilon\_ub}}} \right], \qquad (13) $$

where $M_{\Upsilon\_lb}$ and $M_{\Upsilon\_ub}$ denote the medians of the lower and upper bounds at Ls = $\Upsilon^\circ$, respectively. Meanwhile, $N_{\Upsilon\_lb}$ and $N_{\Upsilon\_ub}$ denote the number of available lower and upper bounds at Ls = $\Upsilon^\circ$, respectively.

### 3. Study Area

Ice blocks analyzed in this study lie at the bottom of an equator-facing steep scarp. The scarp is centered at (85.0° N, 151.5°E) with a total length of ~20 km (refer to its location in Figure 4). This scarp is the same as that studied in Su, Fanara, Zhang, et al. (2023) and termed Scarp 1. The scarp is visually distinguishable by its bright water ice layers (Figure 5). The relatively flat top of the NPLD features a residual layer of mostly dust-free water ice and
can grow to a few meters at most (Figure 4). Over the surface a homogeneous pitted texture at length scales of 10–20 m is revealed, which appears to be the result of differential sublimation (Milkovich et al., 2012; Milkovich & Head III, 2006; Wilcoski & Hayne, 2020). Fallen ice blocks from the steep scarp reside at the surface of the Basal Unit. A typical example of a clustering of ice blocks is shown in Figure 5. The Basal Unit is a low-albedo, interbedded sandy and icy deposits that lie stratigraphically below the NPLD (T. C. Brothers et al., 2015; Grima et al., 2009), and may be one of the largest reservoirs of water-ice on Mars after the PLD (Ojha et al., 2019). The Basal Unit outcrops beneath Scarp 1 belong to the Cavi unit where it transitions into the NPLD (S. C. Brothers et al., 2018; Nerozzi et al., 2022). Although the majority of the ice blocks originate from steep scarps of the NPLD, some of them can also source from the upper part of the Basal Unit with a slope of up to 30° (Su, Fanara, Zhang, et al., 2023).

# 4. Results

## 4.1. Bundle Adjustment, DEM Generation, and Image Orthorectification

A total of 50 images are available for Ice Blocks 1, 2, and 3 with a pixel size ranging from 0.25 m (39) to 0.50 m (11). Most of them (33) were acquired from 2010 to 2013 (MY30–32) with the earliest coverage extended to 2008 (MY29) and the latest ones taken during 2021 (MY36). These images are bundled-adjusted all together. Unfortunately, the bundle adjustment fails to rectify eight images acquired in late winter or early spring when the seasonal layer is thick, due to a lack of small-scale features and “distorted” feature shapes owing to distinct illumination conditions with respect to the rest of the images. Then, the stereo pair ESP_018905_2650 RED and ESP_019222_2650 RED acquired on 8 August 2010 and 2 September 2010 in the summer of MY30 is used to create a DEM (the same pair as used in Su, Fanara, Zhang, et al. (2023)). The stereo pair respectively imaged the Scarp 1 area with spacecraft rolls of −3.6° and −15.5° to achieve a good base-to-height ratio for the purpose of obtaining a reliable 3D model of the surface. The emission angles for these two images are 4.0° and 17.0°, respectively. These images were acquired relatively close in time, with an interval of 25 Earth days, which can minimize differences in solar illumination and surface properties between the stereo pairs. Similarities between the stereo pair can aid in the image matching and parallax determination processes to reconstruct the object heights. The DEM is gridded with 1 × 1 m pixels (Figure 5). Limited interpolation spikes and long lines due to CCD seams have been spotted while inspecting the hillshade generated from the DEM, verifying the correctness and reliability of the bundle adjustment and image matching processes. The bundled-adjusted images are then automatically orthorectified using the aforementioned DEM and used for measurements of the shadows. Visual examination of the orthorectified images show they mutually align well at their overlapping portions. In combination with solar conditions at the acquisition times of the images, local slope and aspect are extracted from this DEM for converting the shadow length variations to temporal thickness evolution of the seasonal CO₂ layer.

## 4.2. Snow/Ice Thickness in Middle-To-Late Spring From THV

We identify three reference ice blocks at Scarp 1 for the purpose of measuring the seasonal snow/ice depth evolution at these specific spots. All of these three ice blocks meet the eligibility criteria. Meanwhile, they are capable of yielding the largest number of reliable and self-consistent thickness measurements among 10 candidates. The other ice blocks are not included in our analysis due to large oscillations of the measured thicknesses there and wider temporal gaps between them. We term the selected three ice blocks as Ice Block 1 (84.984°N, 151.907°E), Ice Block 2 (85.006°N, 151.231°E), and Ice Block 3 (85.008°N, 151.072°E), respectively. Refer to Figure 5 for their locations. The distance between Ice Blocks 1 and 2 is 3,745 m, 832 m between Ice Blocks 2 and 3, and 4,564 m between Ice Blocks 1 and 3. Regional slope at Ice Blocks 1, 2, and 3 are around 11°, 2.5°, and 6°, which are generally gentle compared to that right at the foot of the scarp which can be up to 30°. This could be due to the case that the ice blocks initially falling onto the relatively steep Basal Unit at the foot of the scarps normally can roll or slide further down. Typical HiRISE images of these ice blocks from late winter and spring, when the surface is covered with seasonal deposits, to summer, when the seasonal layer completely sublimates away, are shown in Figure 6. Ice Block 3 is triangular-prism-shaped while the other two feature clear peaks. Ice Blocks 1, 2, and 3 are relatively stable over MY29 to MY36, with average summer heights of ~1.2 m, ~1.4 m, and ~1.8 m, respectively (Section 5.1). During late winter, with the Sun being less than ~5° above the horizon, the CO₂ snow/ice is the thickest. Significant portion of the ice blocks are submerged by the seasonal deposits, and there exists no clear and unique identification of their cast shadows. This situation does not improve much by the very early spring, and only by early to-mid spring, when the deposits gradually sublimate back to the atmosphere, can we
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Upper left: Seasonal snow/ice thickness from MY30 to MY36 for Ice Block 1 at Scarp 1 (85.053°N, 151.833°E) during the northern spring. Bottom left: Seasonal snow/ice thickness evolution from MY29 to MY36 for Ice Block 2. Bottom right: Seasonal snow/ice thickness evolution from MY29 to MY36 for Ice Block 3. Measurements from this study are marked by crosses, and previous Mars Orbiter Laser Altimeter (MOLA) results are marked by dots. It should be noted that for D. E. Smith et al. (2001), the result shown is along the latitudinal annulus centered at 85.5°N. For Aharonson et al. (2004), it is the latitudinal annulus centered at 86°N.</figcaption></figure>

unambiguously distinguish the shadows of the target ice blocks from the background. Continuing into the summer, the surface is entirely free of the seasonal deposits and the surrounding smaller ice blocks reappear.

At Ice Block 1, the images used, solar and slope conditions, along with the ice block heights and seasonal snow/ice depth are summarized in Table S1 in Supporting Information S1. The number of valid measurements during northern spring can vary from three in MY32 and MY36 and five in MY31 to six in MY30. All measurements except for one in MY31 fall within 35° and 70° in solar longitude. Some observations during the springtime have been excluded in the analysis due to bad image quality or no clear identification of the shadows. The thickness evolution patterns of the seasonal deposits at the Ice Block 1 site are shown in Figure 7. These curves generally feature decreasing trends, consistent with the phenomenon that during spring the sun rises above the horizon and the insolation increases with time. It is interesting to note that MY30/32 and MY31/36 share similar patterns, respectively. However, both patterns reach a thickness of ~0.35 m at around Ls = 37.5°. The maximum difference between these two patterns is limited to be within ~0.2 m. The uncertainty of the thickness inversion is about 0.11 m as illustrated in Section 5.1, these interannual differences could thus not be confirmed from a statistical point of view. For Ice Blocks 2 and 3, the images used, solar and slope conditions, along with the ice block heights and snow/ice depth are summarized in Tables S2 and S3 in Supporting Information S1, respectively. Their thickness evolution results are shown in Figure 7. They share similar trends with that of Ice Block 1, although the interannual dispersions are much more obvious than that of Ice Block 1 (but still limited to be less than ~0.2 m). The depth inversion uncertainties at these two ice blocks are 0.10 and 0.16 m, respectively (Section 5.1). Thus, these multiyear variations could not be confidently confirmed. For all of the three ice blocks, the earliest measurements fall within 20° and 25° in solar longitude. There are images acquired in late winter and early spring, but as the elevation angle of the solar rays is less than 10°, the images are dimly lit and attempts to measure the
shadow length have failed. Furthermore, the bundle adjustment of these images can fail and no orthorectified versions are available for the measurements. These issues can be circumvented by the THB approach for which the thickness measurements are extended to cover late winter and early spring (Section 4.3).

Three previous MOLA results in the springtime of MY25, that is, D. E. Smith et al. (2001) at the latitudinal annulus 85.5°N, Aharonson et al. (2004) at the latitudinal annulus 86°N, and Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2022) at the grid element centered at (85°N, 155°E) and of size 0.5° in latitude and 10° in longitude are shown for comparison (Figure 7). Before Ls = 70°, the results from D. E. Smith et al. (2001) and Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2022) are generally consistent with each other, with a maximum thickness of ~1.1 m at the beginning of spring and declining with time. Interestingly, all of the MOLA results consistently predicted a depth of ~0.64 m at around Ls = 70°, before which the thickness measurements of Aharonson et al. (2004) are consistently the largest. When approaching summer solstice (Ls = 90°), all of the MOLA results feature non-zero thickness values, with that of D. E. Smith et al. (2001) being as high as ~0.6 m. These deviations are considered outliers due to biases in the MOLA results as discussed in Section 5.5. Surprisingly, these MOLA results are consistently higher, by a magnitude of up to ~1 m, than our measurements using shadows of the ice blocks. However, it should be noted that at this point the effects of direct condensation onto the tops of the ice blocks have not been considered, and the thickness measurements can be viewed as that contributed solely by snowfalls. These aspects are thoroughly discussed in Section 5.4. Possible explanations of the remaining offsets after the correction for these effects are examined in Section 5.5.

## 4.3. Snow/Ice Thickness in Late Winter and Spring From THB

At Scarp 1, HiRISE images taken in late winter are available only in MY31. Thus, we apply the THB approach to all the images in MY31 with an attempt to place bounds on the snow/ice thickness all the way from late winter to spring, at solar longitudes of 7.0°, 17.3°, 22.7°, 30.3°, 35.4°, 42.8°, 55.1°, 69.5°, and 350.7°, respectively. The corresponding images are ESP_024232_2650 RED, ESP_024509_2650 RED, ESP_024654_2650 RED, ESP_024865_2650 RED, ESP_025010_2650 RED, ESP_025221_2650 RED, ESP_025577_2650 RED, ESP_025999_2650 RED, and ESP_032632_2650 RED, respectively. The solar elevation angles during acquisitions of these images are 7.9°, 12.2°, 14.3°, 17.3°, 19.2°, 21.7°, 25.3°, 28.3°, and 0.9° (late winter), respectively. The results obtained during late winter and early spring serve the purpose of filling the temporal gap of the THV-derived results. Meanwhile, the results obtained during mid-to-late spring can cross-validate the THV-derived results and check if the thickness measurements obtained at Scarp 1 are indeed much lower than the MOLA-derived values. For late winter and early spring images acquired at Ls = 7.0°, Ls = 17.3°, and Ls = 350.7° that fail to be rectified during the bundle adjustment and subsequent orthorectification, we georeference them to the summer image taken at Ls = 128.3° in MY31 (ESP_027674_2650 RED). For this purpose, we select a total of around 10 control points that are homogeneously distributed over the image. These tie points are chosen to be stable ice blocks, corner points of protrusions from layers of the icy stratigraphy in the scarp, and large-scale topographic features over the Basal Unit outcrops. The affine transformation is then applied to tie these image to the summer one when is free of the seasonal layer. We proceed to select the locally tallest ice blocks, that is, with the longest shadows, and locally shortest ice blocks over Scarp 1. After that, the heights of these bounding ice blocks are measured in the summer image and serve as constraints on the seasonal snow/ice depth.

For demonstration purpose, we show example ice blocks used to bound the thickness of the seasonal cover at Ls = 7.0° in early spring and Ls = 35.4° in middle spring in MY31, respectively. At Ls = 7.0°, we display three example ice blocks that have been fully submerged, termed Ice Blocks lb1, lb2, and lb3, which can then indicate the minimum snow/ice depth in the Scarp 1 region (left panel in Figure 8). During the early spring, images reveal no texture related to the underlying ice blocks, demonstrating that the ice blocks have been completely buried underneath. Their heights measured in the subsequent summer when surface is free of seasonal deposits stand at 0.97 m, 1.02 m, and 0.93 m, respectively. After the correction of $-\xi_Y = -0.20$ m, these values become 0.77 m, 0.82 m, and 0.73 m, respectively. We proceed to showcase three other ice blocks that have not been fully submerged, named Ice Blocks ub1, ub2, and ub3, and hence can put upper limits on the snow/ice depth (right panel in Figure 8). Their heights measured at the coming summer stand at 1.20 m, 1.12 m, and 1.10 m, respectively. After the correction of $-\xi_Y = -0.20$ m, these values become 1.00 m, 0.92 m, and 0.90 m, respectively. These lower and upper bounds are extremely self-consistent with mean values of 0.77 and 0.94 m, respectively. Thus, thickness of the seasonal layer at the location of these example bounding ice blocks at Ls = 7.0° is most likely to be within 0.77–0.94 m. At Ls = 35.4° during middle spring in MY31, we show three example ice blocks that have been fully
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Images of example bounding ice blocks (centered at the sub-frames) at Ls = 7.0° in early spring and Ls = 128.3° in summer in MY31, respectively. The images used are ESP_024232_2650_RED and ESP_027674_2650_RED, respectively. These ice blocks are utilized to place bounds on the thickness of the seasonal deposits at Ls = 7.0°. The subscript “lb” means the ice block has been totally submerged and is capable of placing a lower bound on the thickness. Similarly, the subscript “ub” indicates the ice block in question has not been completely submerged and is capable of placing an upper bound on the thickness. Map projection is as in Figure 6. The scale bar applies to all of the images. Illumination is to the bottom-left corner.</figcaption></figure>

submerged, termed Ice Blocks lb4, lb5, and lb6, respectively. Their heights can place lower limits on the snow/ice depth in the Scarp 1 region (left panel in Figure S3 in Supporting Information S1). Their heights measured at the coming summer stand at 0.53 m, 0.61 m, and 0.56 m, respectively. After the correction of $-\xi_Y = -0.25$ m, these values become 0.28 m, 0.36 m, and 0.31 m, respectively. We then illustrate three other example ice blocks that have not been fully submerged, named Ice Blocks ub4, ub5, and ub6, and hence can render indication as to the maximum snow/ice depth (right panel in Figure S3 in Supporting Information S1). Their heights measured at the coming summer are 0.64 m, 0.64 m, and 0.70 m, respectively. After the correction of $-\xi_Y = -0.25$ m, these values become 0.39 m, 0.39 m, and 0.45 m, respectively. These bounding values are largely self-consistent with deviations of less than 0.1 m. Thickness of the seasonal layer at the location of these example bounding ice blocks at Ls = 35.4° is most likely to be within 0.32–0.41 m.

The bounding ice blocks used are quasi-uniformly distributed along Scarp 1 between 150.82°E and 151.95°E in longitude (Figure S4 in Supporting Information S1). Majority of them horizontally fall within 100–600 m away from the foot of the scarp. The number of bounding ice blocks used at each solar longitude ranges from 23 at Ls = 30.3° to 82 at Ls = 17.3° (marked in Figure 9). The resulting correction of $-\xi_Y$ ranges from −0.02 m at Ls = 350.7° to −0.38 m at Ls = 69.5°. The corrected results from THB in MY31, in comparison to that from THV in MY31 and MOLA dataset in MY25, are shown in Figure 9. Note that the THV yield results for Ice Blocks 1, 2,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Results from both THV and THB in MY31, and their comparison to the existing Mars Orbiter Laser Altimeter (MOLA) results. The THV results for Ice Blocks 1, 2, and 3 are marked by crosses and connected by broken lines. The THB results have been corrected for the offset of $-\xi_Y$, and are represented by the violin plots with medians of the upper and lower limits shown as horizontal lines. Upper bounds shown in cyan while lower bounds in orange, and the corresponding number of bounding ice blocks adopted are marked alongside the violin plots. Black squares and dots denote adjusted bounds in consideration of the standard errors of the median constraints. For simplicity and temporal continuity, the constraints obtained at $Ls = 350.7^\circ$ in late winter are plotted at $Ls = -9.3^\circ$.</figcaption></figure>

and 3 (Section 4.2), while the THB results represent large-scale average constraints measured at tens of bounding ice blocks distributed over Scarp 1 (Figure S4 in Supporting Information S1). Ice Blocks 1, 2, and 3, at more than 1.2 m, are much higher than that of the bounding ice blocks, which are at a similar level with depth of the seasonal deposits to be constrained. The THB approach can yield estimates in late winter and early spring while THV are temporally limited to mid-to-late spring. Thickness bounded by median limits are highest during late winter with values most likely falling between 0.92 and 1.02 m at $Ls = 350.7^\circ$, then that decrease to be between 0.63 and 0.65 m at $Ls = 7.0^\circ$ in early spring, between 0.20 and 0.22 m at $Ls = 42.8^\circ$ in middle spring, and gradually decrease to be within 0.01 and 0.03 m at $Ls = 69.5^\circ$ in late spring. The thin layer at $Ls = 69.5^\circ$ is not in conflict with the statement of Piqueux et al. (2015) that the area of the SNPC shrunk to zero at around $Ls = 80^\circ$ in MY31. Differences between upper and lower limits are smaller than 0.1 m at all examined solar longitudes, demonstrating the high precision of the average thickness estimates from THB approach. Apart from shadow length measurement errors (well less than 0.7 pixel which translates to less than 0.1 m in ice block height and thickness bound for image resolution of 0.25 m), errors due to varying solar condition, and the inclusion of overly tall (in the uncovered case) or short (in the completely covered case) bounding ice blocks, dispersion of the thickness constraints at each examined solar longitude can, to some extent, reflect regional variability of the snow/ice depth. The bulk upper and lower constraints overlap well at all examined solar longitudes. This can be partially attributed to, other than the aforementioned factors, the tendency that ice blocks used to derive lower bounds are more likely to be located in areas where snow/ice cover thickens, while that for generating the upper bounds more frequently fall within regions with shallower snow/ice deposits. The adjusted bounds taking into consideration the standard errors of the medians from Equation 13 are represented by black squares and dots in Figure 9, respectively. The updated bounding intervals feature half ranges from 0.03 m at $Ls = 17.3^\circ$ to 0.13 m at $Ls = 350.7^\circ$, and with typical values of ~0.05 m (Table S4 in Supporting Information S1). The measurements from THB are largely consistent with that from THV, especially when considering the uncertainty of the latter which ranges from 0.10 to 0.16 m (Section 5.1). This demonstrates the feasibility and correct implementation of the proposed approaches. Interestingly, both of these results offset the MOLA ones by up to 1 m. However, it should be noted the overlying layers over the ice blocks due to direct condensation have not been considered in the THV and THB models. These thickness measurements thus represent the sole contribution by the snowfalls. Estimation of these effects
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Ice block heights (dots) measured at summertime when free of snow/ice cover and the resultant reference heights (horizontal lines), the latter of which are used for the calculation of the CO₂ snow/ice thickness in the springtime. Note that there exists no height measurements in MY33, MY35, and MY36 due to the lack of image coverage, and the measurement at the adjacent MY34 is used instead as the reference height. Temporal separation between contiguous ticks in the horizontal axis is 90° in solar longitude, that is, a Martian season. su. is the abbreviation for summer.</figcaption></figure>

and the correspondent correction of the THB thickness constraints in MY31 are discussed and described in Section 5.4.

# 5. Discussion

## 5.1. Long-Term Stability of the Ice Block Heights and Precision of the THB Thickness Measurements

Here, we discuss the stability of the heights of Ice Blocks 1, 2, and 3 which are used as reference when computing the snow/ice thickness during the springtime. We also discuss the precision of the resultant thickness measurements from the THB method.

Ice block heights measured when the CO₂ snow/ice cover completely sublimated in various summers are shown as dots in Figure 10. Summer images used are marked in bold in Tables S1, S2, and S3 in Supporting Information S1 for Ice Blocks 1, 2, and 3, respectively. When multiple ice block height measurements exist within the same summer, their average is taken as the reference value (horizontal lines in Figure 10). The bare heights of Ice Blocks 1, 2, and 3 are generally stable with variations within ~0.1 m for Ice Blocks 2 and 3 but can be up to 0.2 m for Ice Block 1. It is worthy to note that these heights generally exhibit a gentle decline from MY29 to MY33, after which an increase in MY34 can be spotted. These long-term decreases in ice blocks heights could be due to tumbling under its own gravity or that triggered by external forces. Meanwhile, aeolian erosion and ice sublimation may also be responsible. To mitigate the effects of quasi-stability, we assign the heights measured at temporally adjacent summers as the baseline for inverting for thickness of the seasonal deposits during the spring (Figure 10). While it is reasonably to speculate their decreases in heights due to wind erosion and disintegration of the upper parts due to condensation/sublimation of the seasonal ice cover, we hardly find any physical explanations for their growth. Thus, we attribute these abnormal height increases to potential errors in the reference measurements. The induced systematic biases in the thickness values obtained in the springtime for Ice Blocks 1, 2, and 3 are then 0.06 m, 0.05 m, and 0.07 m, respectively.

To look into the random errors that can also deteriorate the precision of our thickness measurements, we combine two types of observables. If there are two or more height measurements during each Martian summer, the inconsistencies between measurement pairs are used as a reasonable gauge of the random errors. Additionally, we examine the repeatability between thickness measurement pairs in the springtime if their temporal separations are within an empirical threshold of 5° in solar longitude. In this case, seasonal snow/ice sublimation within these
relatively short timespans can be neglected and the aforementioned deviations from zero can gain us insights into the precision of our measurements. The standard deviation is calculated as follows:

$$ \sigma_H = \sqrt{\frac{1}{n-1} \sum_{i=1}^{n} \left(\frac{\Delta H_i}{2}\right)^2}, \qquad (14) $$

where $\Delta H_i$ represents the ith qualified height discrepancy, out of a total of n, to be evaluated. We scale the random error $\sigma_H$ by a factor of three for the metric to be robust. For Ice Block 1, there are three qualified values during the summers and seven values during the springs. The random error of Ice Block 1 is 0.09 m. For Ice Block 2, there are five qualified height discrepancies during the summers and nine values during the springs, respectively. And, the random error of Ice Block 2 is calculated at 0.09 m, which is the same as that of Ice Block 1. Meanwhile, for Ice Block 3, five qualified values are available during the summertime and nine during the springtime. The random error of Ice Block 3 stands at 0.15 m.

We proceed to calculate the Mean Square Error (MSE) that takes into consideration of both the systematic bias expressed as $\epsilon_H$ and random error, expressed as $3\sigma_H$, at each of the three locations. These quantities are related as

$$ \text{MSE}(\hat{\Omega}) = \mathbb{E}\left[\left(\hat{\Omega} - \Omega\right)^2\right] = \text{Var}(\hat{\Omega}) + \text{Bias}^2(\hat{\Omega}) = 9\sigma_H^2 + \epsilon_H^2, \qquad (15) $$

where $\Omega$ denotes the thickness measurement during springtime. The square roots of the MSE errors for thickness measurements for Ice Blocks 1, 2, and 3 are 0.11 m, 0.10 m, and 0.16 m, respectively. The uncertainty of our thickness measurements from ice blocks can thus be ~0.1 m at the best scenarios, which is on similar level with the statement from Mount and Titus (2015). These empirical uncertainties represent the cumulative impacts of a combination of factors such as variations in image quality, uncertainties in HiRISE DEMs, image orthorectification errors, and shadow length measurement errors. We note that changing solar elevation and azimuth can also possibly lead to varying cross-sections of the ice blocks perpendicular to the solar illumination and introduce errors into the shadow length measurements. As Ice Blocks 1 and 2 are lone-peaked, their illuminated tips will remain largely unchanged under varying solar condition. Instead this factor could play a role for Ice Block 3 which is triangular-prism-shaped.

Although these precision statistics do not allow unambiguous confirmation of possible seasonal thickness variations over multiple years (up to ~0.2 m), they do validate the notion that no significant thickness variations exist in middle-to-late springs from MY29 to MY36. In addition, it strengthens the statement that our thickness measurements are much lower than previous MOLA results. This statement is independently verified using the THB approach as shown in Section 4.3.

## 5.2. Characteristics of the THB Approach

Compared to THV, the THB approach avoids the need to clearly identify and accurately measure the shadow length of the ice blocks during late winter and spring. Instead, one only needs to check if there exist cast shadows around the ice blocks to decide whether they have been completely covered or not. That is a huge advantage over the THV approach, especially during late winter or early spring when the shadow boundaries of the ice blocks can be unidentifiable due to low solar elevation and thick seasonal layer. This advantage means that images obtained during wintertime and springtime are not necessarily required to be bundled adjusted and orthorectified, which can greatly boost computation efficiency. In fact, bundle adjustment can fail with images acquired in late winter or early spring which cannot be utilized in the THV approach. For these aforementioned reasons, the THB approach can temporally extend to cover late winter and early spring. Unfortunately, THB can be difficult to implement during late spring as the snow/ice depth shrinks to a low level. In such case, the bounding ice blocks along with their shadows required to place the lower and upper limits become too small to be resolved in the HiRISE images. In addition, the strength of the constraints that can be placed by the THB approach depends on the height gap between the bounding ice blocks. For scarps with good concentration of ice blocks, the precision of the obtained thickness constraints would be satisfying. Further improvement can be expected by spatial averaging of the thickness bounds over a large area. However, THB can fail in regions where number of ice blocks are limited or the distribution of the ice block heights does not overlap with the seasonal snow/ice thickness to be constrained,
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison of thickness results in early springs of MY31 (same as in Figure 9), MY32, and MY36 from the THB approach. The bounding limits have been corrected for the offset of $-\xi_Y$. The results are represented by the violin plots with medians of the upper and lower limits shown as horizontal lines. The corresponding number of bounding ice blocks adopted are marked alongside the violin plots. Black squares and dots denote adjusted bounds in consideration of the standard errors of the median constraints (Equation 13). For simplicity and temporal continuity, the constraints obtained at $L_s = 350.7^\circ$ in late winter of MY31 are plotted at $L_s = -9.3^\circ$. Existing Mars Orbiter Laser Altimeter (MOLA) results are shown for reference.</figcaption></figure>

for example, in fields where only debris and extremely small ice blocks remain. Despite various merits of THB over THV, these two approaches can be complementary to each other. For temporal extension to late winter and early spring, and spatial extension to a large set of ice block clusterings at individual scarps, automatic version of THB would be the choice (Section 5.8). However, in case high spatial resolution is required within regions of interest, for example, over specific surface features, or when available ice blocks are scarce, THV should be applied instead. Meanwhile, in a situation when both approaches can work, their results at mid-to-late springtime can be used for cross-validation purposes.

### 5.3. Interannual Snowfall Variations From THB

Due to high precision of the THB method over scarp-scale area, we proceed to apply it to examine possible interannual thickness variations in the very early spring (before $L_s = 10^\circ$) when the seasonal layer is at its thickest. After excluding images with bad quality, only ESP_024232_2650_RED acquired at $L_s = 7.0^\circ$ in MY31, ESP_033120_2650_RED at $L_s = 9.7^\circ$ in MY32, and ESP_068224_2650_RED at $L_s = 3.5^\circ$ in MY36 are available. For MY31 and MY32, the same summer reference image is used (ESP_027674_2650_RED) and the corrections of $-\xi_Y = -0.20$ m and $-\xi_Y = -0.22$ m are respectively applied to the overestimated thickness values. For MY36, the summer image acquired in MY34 is used (ESP_053730_2650_RED) and a correction of $-\xi_Y = -0.08$ m is adopted. Bounding ice blocks are selected over the entire Scarp 1. The corrected results and their comparison are shown in Figure 11. The results at $L_s = 350.7^\circ$ and $L_s = 17.3^\circ$ in MY31 are also shown alongside that at $L_s = 7.0^\circ$ in MY31 to illustrate the general depth evolution trend during the beginning of spring (same as in Figure 9). The median depth bound at $L_s = 9.7^\circ$ in MY32 is from 0.74 to 0.76 m. Meanwhile, the median depth bound at $L_s = 3.5^\circ$ in MY36 ranges from 0.98 to 1.03 m, which is ~0.36 m larger than between 0.63 and 0.65 m at $L_s = 7.0^\circ$ in MY31. Considering the adjusted bounds at the examined solar longitudes in MY31 and MY36 (black squares and dots in Figure 11), the interannual variation of ~0.36 m can be deemed credible. Indeed, treating the thickness bounds at MY31 and MY36 as being independent, the propagated uncertainty associated with the MY31 to MY36 thickness variation is 0.13 m. A supporting evidence is that there exist plenty of ice
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Images of example ice blocks (centered at the sub-frames) used to deduce depth of topping crowns. Acquisition time stamps of the images are Ls = 42.8° in middle spring and Ls = 128.3° in summer in MY31, respectively. The images used are ESP_025221_2650_RED and ESP_027674_2650_RED, respectively. Left panel shows an ice block with parallel walls and its shadow, the thinning of which from middle spring to summer is obvious. Right panel refers to the change of a cone-shaped ice block. Map projection is as in Figures 6 and 8, and Figure S3 in Supporting Information S1. The scale bar applies to all of the images. Illumination is to the bottom-left corner.</figcaption></figure>

blocks that are visually detectable at Ls = 7.0° in MY31 can stay completely covered at Ls = 3.5° in MY36. We can also see in the plot that all of the examined thickness values are lower than the existing MOLA results, by a magnitude of 0.1–1 m. However, we note again that the effects of crowning layers on top of the ice blocks due to direct condensation have not been considered in the applied THB model. As such, the detected multiyear thickness variation refers to that solely induced by the snowfalls (Section 5.4). This is the first time that an interannaul variation in the amount of Martian snowfalls is detected, stressing the importance of carrying out long-term monitoring of the Martian volatile and dust cycles.

## 5.4. Crown Depth Estimation and Correction in MY31

In the presented geometric models (Figures 1–3), we have not considered the effects of moating and crowning. For rocks, moats form in the fall while they are still warm with stored heat energy from summertime insolation (Mount & Titus, 2015). Snow/ice falling around the warm rocks will sublime, leaving a void. Moating has been extensively observed around the rocks in the dune fields (Mount & Titus, 2015), which can be attributed to the large contrast in thermal inertia of rocks (relatively high) and sands (relatively low). The ice blocks fallen from the NPLD, as investigated in our case, are made mainly of water ice, which feature very high thermal inertia. However, the upper layer of the Basal Unit deposits where the ice blocks reside is presumably composed of previously fallen ice blocks, leading to its sloped surface at the foot of the NPLD. Evidence has it that the potentially active scarps generally feature a higher overall slope from Basal Unit outcrop top to bottom than the potentially quiescent ones (Russell et al., 2010). As a result, the contrast in thermal inertia between the ice blocks and the underlying Basal Unit top layer should be limited. Thus, extensive moating well into the springtime should be unlikely. Indeed, we have not observed any instance of moating around the ice blocks during springtime with typical examples shown in Figures 6 and 8, and Figure S3 in Supporting Information S1. However, moating in late northern spring due to radiation of the ice blocks can be possible, which will be discussed in Section 5.6. The formation of a topping ice layer, that is, crowning, has been properly accounted for in Mount and Titus (2015) when using shadowing of the rocks to invert for thickness evolution of the seasonal snow/ice deposits. For rocks, crowns form in the winter, after they have lost sufficient heat and are expected to be most prevalent in early to mid-spring. If crowns form then the measured shadow length will be longer than in the case without, making the uncorrected rock height larger and thus reducing the calculated snow/ice thickness. For ice blocks, the crowning can also happen and should be properly accounted for. Indeed, we have observed systematic widening of the ice block walls in the HiRISE images which is due to direct condensation of CO₂ as frosts (Figure 12).

Mount and Titus (2015) corrected for the crowning and moating by adding the increases in snow/ice depths between consecutive solar longitudes to all prior depths. This forces the ice depth to either decrease or remain the same between different solar longitudes, and the corrected values serve as minimum bounds to the seasonal snow/
ice depth. The representative ice depth curve before correction in their study features significant increases during middle-to-late spring (Figure S3 in Mount & Titus, 2015). However, in our case we have not seen significant off-season increases (Figures 7, 9 and 11), so the results should not change much after incorporating their correction scheme. This indicates that the depths of the crowns over the ice blocks in this study likely quasi-linearly decrease to zero right at when the seasonal deposits completely sublimate back into the atmosphere. To conclude, the correction scheme described in Mount and Titus (2015) can be inapplicable in our case. We thus propose an original way to estimate and correct for the effects of the crowning on top of the ice blocks. The proposed correction scheme relies on three reasonable assumptions: (a) Atmospheric deposition as snowfalls do not accumulate over the top of the selected spiky ice blocks and over steep ice block walls. In fact, it is known that aeolian processes can redistribute the carbon dioxide snow crystals into topographic lows, damping the topography (Mount & Titus, 2015). Ice blocks shaped like triangular prisms or those that feature lone peaks are used in the demonstration of the THV and THB approaches (Figures 6 and 8, Figure S3 in Supporting Information S1), which makes significant accumulation of snow particles over the ice blocks unlikely. Despite the thin Martian atmosphere (roughly 1% thinner than that on Earth), winds on Mars are more than capable of moving sand and fine particles of snow. Moreover, strong katabatic winds can be expected at the foot of these steep scarps. Thus, under this assumption the crowns in our case can be solely composed of directly condensed frosts; (b) Frosts through direct condensation are homogeneously distributed around the ice blocks and their adjacent surroundings. On Mars, CO₂ constitutes 95% of the total atmospheric pressure and condensates mainly as compact slab ice with low porosity (Grisolle, 2013; Portyankina et al., 2019). A good terrestrial analogy is the glaze ice or glaze frost, which is a smooth, transparent, and homogeneous ice coating that forms when freezing rain droplets hit a surface. This hypothesis can be corroborated by the fact that there exist limited thermal contrast between the ice blocks and the Basal Unit outcrops where they reside. It is further justified by similar depths of ice block crowns at a specific solar longitude over the entire scarp, as illustrated in Figure 13; (c) Negligible moating around the ice blocks, as previously demonstrated in this section. These assumptions mean that the contributions of the snowfalls and frosts to the thickness of the seasonal ice deposits can be decomposed and estimated separately. The thickness measurements from THV and THB (ΩY at Ls = γ°) can be considered as merely contributed by snowfalls (Figures 7, 9 and 11). Indeed, if there exists no snowfall, then the shadow length of the ice blocks would remain unchanged and thickness measurements from THV and THB should all equal zero. These assumptions also enable us to use the widening of the ice blocks as a proxy to approximate the depths of the crowning frost layers over the ice blocks (ΩY at Ls = γ°). The total thickness of the seasonal layer then can be calculated as the sum of ΩY and ΩY. The geometry of the ice blocks and seasonal deposits after the assumptions made is shown in the top panel of Figure S7 in Supporting Information S1. We explore possible violations of these assumptions and their impact on our snowfall thickness measurements in Section 5.6.

To approximate the depths of the crowns, we first examine the widening of cylinder-shaped ice blocks with parallel walls, mostly vertical walls, during late winter and spring as compared to their dimensions during the summertime when the enveloping seasonal CO₂ frost layer completely sublimates away. A typical example of these particular ice blocks is shown in Figure 12. Within the portion of the shadow with parallel bounding lines, the width of the shadow, and hence that of the ice block itself, can be directly measured along any transverse line running perpendicular to the solar illumination. Then, the crown depth at Ls = γ° can be approximated by

$$ \Omega_Y^c = \frac{W_Y^1 - W^2}{2}, \qquad (16) $$

where W<sub>Y</sub><sup>1</sup> is the width of the shadow of the ice block during late winter or spring, and W<sup>2</sup> denotes that of the same ice block during summertime. Temporal evolution of the crown depth over the late winter and spring in MY31 is illustrated in Figure 13. The abstract probability distribution of crown depth measurements at a specific solar longitude is approximated using Equation 11 and represented using a pink violin. At Ls = 42.8°, the image features the best sharpness and the most measurements are made (a total of 23), with a median value of 0.26 m. In contrast, at Ls = 350.7° in late winter the image is the most severely blurred and only one valid measurement is obtained with a value of 0.64 m. For early spring images acquired at Ls = 7.0° and Ls = 17.3°, the dispersions of the measurements are relatively large due to degraded quality associated with still low solar elevation angles. As can be seen in the plot, the crown depth quasi-linearly decreases from 0.64 m at Ls = 350.7° in late winter, 0.26 m at Ls = 42.8° in middle spring, to 0.045 m at Ls = 69.5° in late spring. The approximately quasi-linearly declining trend corroborates our previous conjecture in the light of limited middle-to-late spring thickness increases present
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Depth estimates of the crowns from direct condensation over Scarp 1 in MY31, as represented by the violin plots with medians shown as horizontal lines. The corresponding number of independent measurements are marked alongside the violin plots. Pink violins represent measurements using ice blocks with parallel walls and the blue one at Ls = 42.8° refers to that from ice blocks with nonparallel walls. Black and blue dots denote derived bounding intervals by taking into account of the standard errors of the medians. For simplicity and temporal continuity, the single constraint obtained at Ls = 350.7° in late winter is plotted at Ls = -9.3°.</figcaption></figure>

in our uncorrected results. Like it has been done in Equation 13, we derive adjusted bounding intervals for the most likely crown depths at Ls = Y° by taking the standard errors of the medians into account:

$$ \Theta_Y^c = \left[ M_Y^c - \frac{3 \times \text{SMAD}_Y^c}{\sqrt{N_Y^c}}, M_Y^c + \frac{3 \times \text{SMAD}_Y^c}{\sqrt{N_Y^c}} \right], \quad (17) $$

which are marked as black dots in Figure 13. The half ranges of these bounding intervals are from 0.02 m at Ls = 55.1° to 0.09 m at Ls = 7.0°, with typical values on the order of 0.05 m (Table S4 in Supporting Information S1). The limited dispersions in crown depths, or ice block widening, over the entire scarp at individual solar longitudes can to some extent reflect the spatial uniformity of the directly condensed layer in the vertical dimension. This serves as an additional line of evidence to justify our second assumption that the directly condensed frost layer is of the same thickness over the ice blocks and their adjacent surroundings. The estimated 0.64 m thick frost layer due to direct condensation in late winter is significantly shallower than the snow layer with a thickness of 0.97 ± 0.13 m (Section 4.3). In terms of thickness and volume, the snowfalls can make up 60.2% of the seasonal deposits in late winter. However, assuming the crowns to be slab ice that features a density greater than or equal to 1,190 kg/m³ (≤26% porosity, ≥8,000 µm grains) and the snow to feature a density of less than or equal to 420 kg/m³ (≥74% porosity, ≤1,000 µm grains) (Mount & Titus, 2015), then snowfalls should account for less than 34.9% of the seasonal deposits in terms of mass during late winter.

As majority of the ice blocks used in the demonstration of the THV and THB methods are peaked or triangular-prism-shaped, we wonder if similar amount of frosts is directly condensed onto these nonparallel walls. Thus, we also look into the widening of the cone-shaped ice blocks. As the boundaries of the shadows of these peaked ice blocks are triangular-shaped, shadow width measurements have to be made at properly adjusted positions so that they can be directly compared to reveal the widening. The geometry and the procedure to calculate unbiased widening of the cone-shaped ice blocks in spring are illustrated in Text S1 in Supporting Information S1. This approach requires clear prior identification of the triangular-shaped shadow boundary and thus permits many fewer valid estimates than those using the ice blocks with parallel walls. In fact, it is inapplicable to images
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Thickness evolution of the seasonal layer in MY31 from THB after the correction for the crowns over the ice blocks, and their comparison to the existing Mars Orbiter Laser Altimeter (MOLA) results. Uncertainty bars denote the accumulated errors from both the uncorrected thickness measurements and the crown depth estimates. For simplicity and temporal continuity, the constraints obtained at Ls = 350.7° in late winter of MY31 are plotted at Ls = -9.3°.</figcaption></figure>

acquired in late winter and early spring. Here, we only apply this approach to the image acquired at Ls = 42.8° which features the best quality in MY31. We obtain a total of 11 valid measurements that range from 0.22 to 0.36 m, with a median of 0.29 m (blue violin in Figure 13). One of these measurements reaches convergence after two iterations while all the rest satisfy the stopping criterion after just one iteration. These measurements overlap well with those obtained using ice blocks with parallel walls with a median of 0.26 m. Meanwhile, the bounding interval in light of the standard error of the median constraint is highly consistent with that from the ice blocks with parallel walls (blue and black dots in Figure 13, respectively). That means the steepness of the ice block walls does not significantly influence the depth of the directly condensed frosts. This conclusion also corroborates our second assumption that the thickness of the directly condensed frosts is the same over the ice block walls and top, which permits us to use the widening of the ice blocks as a proxy to approximate the depths of the overlying crowns.

Eventually, we correct the snowfall thickness measurements from THB in MY31 for these direct condensation effects (Figure 14). The total thickness of the seasonal deposits at Ls = Y° is calculated as follows:

$$ \Omega_{Y}^{\text{sum}} = \Omega_{Y} + \Omega_{Y}^{c} = \frac{\bar{\Theta}_{Y}^{\Omega} + \bar{\Theta}_{Y}^{c}}{2}, \qquad (18) $$

with its bounding interval set by treating the snowfall thickness estimates and the crown depths as two independent variables:

$$ \Theta_{Y}^{\Omega^{\text{sum}}} = \left[ \Omega_{Y}^{\text{sum}} - \frac{\sqrt{\Delta\Theta_{Y}^{\Omega^{2}} + \Delta\Theta_{Y}^{c}}}{2}, \quad \Omega_{Y}^{\text{sum}} + \frac{\sqrt{\Delta\Theta_{Y}^{\Omega^{2}} + \Delta\Theta_{Y}^{c}}}{2} \right], \qquad (19) $$

where $\bar{\Theta}_{Y}^{\Omega}$ and $\Delta\Theta_{Y}^{\Omega}$ denote the average and range of the bounding interval for the snowfall thicknesses from the THB approach (Equation 13 and Figure 9), respectively. Meanwhile, $\bar{\Theta}_{Y}^{c}$ and $\Delta\Theta_{Y}^{c}$ are the average and range of the crown depth bounding interval by examining the ice blocks with parallel walls (Equation 17 and Figure 13), respectively. For Ls = 350.7° in late winter, there exists no bound interval for the crown depth as there is only one
valid estimate. We then tentatively set its interval range to be two times that at the proximate $L_s = 7.0°$. The bounding intervals from Equation 19 are shown as error bars in Figure 14. The corrected thickness stands at 1.63 m with a half interval range of 0.22 m at $L_s = 350.7°$ in late winter, 0.45 m with a half interval range of 0.06 m at $L_s = 42.8°$ in middle spring, and decreases to 0.06 m with a half interval range of 0.05 m at $L_s = 69.6°$ in late spring. A majority of the bounding interval half ranges are less than 0.1 m (Table S4 in Supporting Information S1). The offsets between these corrected thicknesses and existing MOLA results are thus significantly reduced. At late winter, our thickness estimate of $1.63 \pm 0.22$ m is ~0.5 m above that of D. E. Smith et al. (2001) and Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2022). However, it is compatible with that of Aharonson et al. (2004) considering the associated uncertainty of ~0.2 m. Then, the thickness estimates in this study quasi-linearly decline toward summer solstice. In comparison, the MOLA results decline more gently than ours before $L_s \approx 60°$ in late spring. The gap between our results and the MOLA ones thus forms at the beginning of spring and enlarges to be up to 0.8 m at late spring. These lingering offsets and possible reasons behind are thoroughly discussed in the subsequent Section 5.5.

From late winter to spring, the solar elevation angle increases and so does the insolation. If assuming the albedo and average density of the seasonal deposits to be constant throughout winter and spring, then the slope of the thickness evolving curve should increase with time. The relatively constant declining rate of our curve as in Figure 14 thus indicates the combined effects of density variations, like densification due to gravity-induced self-compaction and re-crystallization (Eluszkiewicz et al., 2005; Mount & Titus, 2015; Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022), and albedo variations that take place in the seasonal deposits (Gary-Bicas et al., 2020; Pommerol et al., 2013). Indeed, we note that thickness evolution of the seasonal deposits at the three study sites in Mount and Titus (2015) features distinct trends during the spring: quasi-constant thinning rate at Phoenix's landing site (68°N, 233°E); temporally decelerating thinning rate at a dune field labeled Dunes (75°N, 282°E); temporally accelerating thinning rate at a crater site called Louth (70°N, 103°E).

## 5.5. Comparison to Previous MOLA Results

Inconsistencies up to 0.4 m during the northern spring exist between the existing MOLA results (Figures 7 and 9). The difference between the curve of Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2022) and that of Aharonson et al. (2004) could mainly be attributed to different procedures utilized to post-correct for the global temporal bias in the MOLA dataset. This bias has been observed at all latitudes, and it features a peak-to-peak amplitude of 2 m and a phase that resembles the synodic period of Mars (Xiao, Stark, Steinbrügge, et al., 2022). Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2022) assumed the bias to be constant at the north polar region, obtained it at annulus 50°N, and used it to correct the thickness curves retrieved at the polar regions. In contrast, Aharonson et al. (2004) parameterized the 3D systematic offsets of the MOLA altimetric profiles by overlapping polynomials and solved for the coefficients by minimizing height residuals at cross-overs located between 57°S/N together with that acquired within 15 Earth days in the polar regions. However, the cross-overs in the polar regions used in their adjustment can only help to improve the self-consistency of the profiles there. Thus, their correction for this global temporal bias at the polar regions was essentially extrapolated from the bias obtained within the equatorial regions. As for D. E. Smith et al. (2001), they treated the thickness variations acquired at 60°S/N annuli as the bias and corrected for it at each polar annulus. Their thickness variations are surprisingly lower compared to the other MOLA results, which may arise from over-smoothing they have applied. Unfortunately, Xiao, Stark, Steinbrügge, et al. (2022) confirmed latitude-dependence of the global temporal bias, and by treating it as a constant across the entire polar regions can introduce temporal bias into the MOLA thickness evolution curves. Thus, it is imperative that we locate the physical cause for this bias and accurately model it for the purpose of correcting the MOLA heights.

After the correction for direct condensation effects, offsets up to ~0.8 m in magnitude between the results using shadowing of the ice blocks and the available MOLA results can still be observed during the spring of MY31 (Section 5.4). Here we discuss several aspects that could possibly be responsible for these offsets: (a) Discrepancies in geographical and year-to-year coverage that exist between the results. Result from D. E. Smith et al. (2001) represents the average pattern at latitudinal annulus 85.5°N, and that of Aharonson et al. (2004) at latitudinal annulus 86°N. Meanwhile, result from Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2022) is representative of thickness evolution at a grid element of size 0.5° in latitude and 10° in longitude. Furthermore, MOLA results date back to MY24/25 while results in this study are from MY29 to MY36. Therefore, geographical differences and inter-annual variability might partially contribute to the offsets; (b) MOLA inherent
biases. An image mosaic at the north polar regions of Mars has been made by Wang and Ingersoll (2002) using the Mars Global Surveyor (MGS) Mars Orbiter Camera (MOC) wide-angle swaths at Ls = 111° in MY25. No apparent seasonal deposits can be observed at that time, which is right after the summer solstice, outside of the NPLD. We also search through the MOC dataset and compare the wide-angle swath acquired at Ls = 75.7° to that taken at Ls = 83.2° during the northern spring in MY25 (Figure S5 in Supporting Information S1). We observe that the seasonal deposits which were present at Ls = 75.7° at Scarp 1 had largely sublimated away before Ls = 83.2°. Meanwhile, Piqueux et al. (2015) carried out continuous tracking of the SNPC edges at multiple Mars Years using MGS Thermal Emission Spectrometer and showed that it completely sublimated away at around Ls = 80° in late spring of MY25. Thus, the biases are at least ~0.60 m (D. E. Smith et al., 2001), ~0.35 m (Aharonson et al., 2004), and ~0.41 m (Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022) for these MOLA results, respectively. In fact, various complicated processes can bias the MOLA results, for example, pulse saturation due to high albedo of the seasonal deposits, the lack of a variable gain amplifier, and limited recording digital ranges (Neumann et al., 2001, 2003), interference of the laser pulses with dynamic seasonal features, incomplete correction for the global temporal bias, and penetration of the laser pulses into the translucent slab ice (Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022). Unfortunately, each of these listed factors is difficult to model and quantify; (c) Interannual variations. Possible variations in the quantities of snowfalls and frosts can exist between MY24/25 (acquisition time of MOLA footprints) and MY31 (e.g., Figure 11). However, these multiyear variations should generally decrease with time in spring as the seasonal layer thins and finally vanishes, which contradicts the observed temporal evolution of the offset; (d) Lower snow/ice thickness expected at the ice block fields than that in areas further away from the scarps. The ice blocks used for the purpose of this study are located within 1 km of the steep equator-facing scarps (Figure S4 in Supporting Information S1). Areas closer to these scarps are expected to retain and absorb more heat relative to areas further away. Cooling during the northern fall/winter in the ice block fields could be less efficient due to reduced sky-view and therefore reduced ability to lose heat via infrared radiation to space. Then, during the northern spring/summer increased incident radiation should be expected in the ice block fields due to infrared radiation being emitted from the nearby equator-facing scarps and the steep BU slopes. Even slightly increased temperature can significantly enhance the snow/ice ablation rates. Therefore, a lower seasonal deposit thickness should be expected at the study region than in areas further away from the scarps. To conclude, there exists strong suggestion to argue that a large portion of these discrepancies result from biases inherited in the MOLA dataset.

## 5.6. Comparison to Predicted Snowfall Quantities

In this section, we discuss our measured snowfall thickness as compared to that predicted by a simple CO₂ cloud settling model using atmospheric profiles acquired by MCS onboard the MRO (Alsaeed & Hayne, 2022). We now know that the snowfall rates are significantly larger in the north than those in the south, because the former features larger surface air pressure and a thicker water ice polar hood that can serve as condensation nuclei for the CO₂ snowfalls (Alsaeed & Hayne, 2022; Gary-Bicas et al., 2020). Alsaeed and Hayne (2022) concluded that the CO₂ snowfall equivalent thickness accumulated throughout fall and winter is on the order of several millimeters poleward of 65°N. Surprisingly, these theoretical values are two orders of magnitude smaller than the measured ~0.97 ± 0.13 m thick snow layer during late winter at Scarp 1. The depth of ~0.97 ± 0.13 m due to snowfalls even significantly exceeds that of the frost layer due to direct condensation (0.64 ± 0.18 m at late winter). However, it is worth mentioning that there exist gaps in the MCS data in the lowest atmosphere (approximately up to 5 km in altitude) where the instrument cannot probe through the optically thick clouds. This means that the effects of localized storms which can rapidly drive up snowfall rates were not captured in Hayne et al. (2014); Alsaeed and Hayne (2022). In addition, zonal averaging applied in the MCS analysis can further erase out aggregation of the snowfalls at local scales. Thus, the theoretical thickness values from Alsaeed and Hayne (2022) can be considered as the lower limits. In fact, it is still unknown which depositional mode (snowfall vs. frost) prevails in the north polar region of Mars. A previous snowing model applied to the southern hemisphere indicated that snowfall can contribute between 3% and 20% by mass of the seasonal deposits (Hayne et al., 2012, 2014). Colaprete et al. (2005, 2008) demonstrated that topography-driven atmospheric waves may generate highly heterogeneous but repeatable snowfalls. Thus, how the steep scarps affect the snowfall amount at their foots where the ice blocks lie also needs to be examined. It should be cautioned that in our measures of the snowfall thickness we have assumed no, or negligible, snowfall to be accumulated on top of the spiky ice blocks. However, the mechanical behavior of Martian dry ice snowfall is still largely unknown and its shedding or sluffing over ice blocks can be a complicated process. A variety of other factors, for example, wind speed and properties of the ice blocks
themselves, can also affect the adherence of the snowfall onto them. We have assumed that the directly condensed
frost layer is homogeneous in thickness around the ice blocks. However, different rates of accumulation and
sublimation, aeolian erosion, and metamorphism could possibly lead to crown depth that does not equal that of the
frost layer at the ice block walls. When the solar elevation angle is lower during late winter and early spring than
in late spring, the moats are more likely to cast circumferential shadows around the ice blocks and be observed in
the HiRISE images. Since we have not observed these shadows, moating is unlikely to be pervasive and extensive
during late winter and early spring. As ice blocks receive more insolation during late spring and begin to radiate
heat into the surroundings, forming of moats can be expected. However, these moats should not reach all the way
down to the BU surface otherwise they would be visible in the images. Therefore, moating could probably have an
impact on the derived snowfall thickness during late spring. In Text S2 in Supporting Information S1, we briefly
explore how the violation of these assumptions can bias the snowfall thickness measurements. We show that
unconsidered snowfall remaining on top of the ice blocks can lead to underestimation of the snowfall thickness.
Meanwhile, we show that if the crown is actually thinner than the frost thickness measured at the ice block walls,
then our approaches can overestimate the snowfall thickness. As for the moating during late spring, our ap-
proaches would underestimate the snowfall thickness if the moating is laterally extensive (at least 2 m and ~0.7 m
in radius for THV and THB, respectively). Unfortunately, these biases due to our assumptions could not be
quantified at the current stage. Future simulations of the aforementioned processes within laboratories that can
mimic the Martian climate would be needed. We plan to obtain snowfall thicknesses at active scarps over the
entire north polar region (Section 5.7) and check if our measurements are consistently way higher than the
modeled ones. If so, that would mean the snowfalls on Mars are likely much more frequent and violent than
thought and the cloud settling model applied by Alsaeed and Hayne (2022) may be missing some important
ingredients. As the condensed water ice particles serve as condensation nuclei for the CO₂ particles during at-
mospheric deposition, higher CO₂ snowfall rates then mean a larger fraction of the Martian water cycle happens
through this scavenging mechanism at the poles.

5.7. Availability of Ice Blocks Across the North Polar Region

Here we discuss the spatial availability of the ice blocks as to examine the maximum coverage of upcoming mass-
application of the proposed approaches.

Russell et al. (2010, 2012) examined HiRISE images covering ~70 Basal Unit outcrops around the NPLD for scattered ice blocks and debris on Basal Unit slopes, and qualitatively grouped a total of 20 potentially active outcrops into three categories in terms of “likelihood of recent mass-wasting activity”. None of the peripheral scarps without Basal Unit exposure features likely mass-wasted NPLD detritus at their base (Russell et al., 2012). This indicates an important role the Basal Unit may play in steepness of the scarps and the related mass wasting processes. These active scarps are spatially limited to be within 80°N and 85°N. In addition, there exists no active scarps in Gemina Lingula and Gemini Scopuli due to the lack Basal Unit outcrops (Fishbaugh & Head III, 2005; Nerozzi et al., 2022). However, the number of HiRISE observations available was relatively small back then in 2012.

We seek to carry out a comprehensive search for scarps with ice blocks where our approach can be potentially applicable to get insights into the condensation/sublimation cycle of the polar seasonal deposits. The search is done by merely examining the available HiRISE images. Note that there also exist other optical images captured by Mars Express's High Resolution Stereo Camera (HRSC; Neukum & Jaumann, 2004), the MRO Context Camera (CTX; Malin et al., 2007), and the ExoMars Trace Gas Orbiter (TGO) Colour and Stereo Surface Imaging System (CaSSIS; N. Thomas et al., 2017). However, the best resolution of these images can range from 4 to 12.5 m/pixel, which is not enough to even capture the ice blocks with maximum size of ~5 m (Fanara et al., 2020b). The search combines two aspects of efforts: (a) We locate all of the sites over the NPLD and their close vicinity which have been observed for at least five times, and randomly select a summer image for visual inspection of fallen ice blocks. Multiple images of the same sites can facilitate the cross-validation of the obtained snow/ice thickness values and enable insights into the interannual depth variations (Section 5.3). It should be noted that these regions of interest are not limited to steep marginal scarps, but can also cover spiral troughs that expose stratigraphic layers (Figure 4), ice-filled craters over the residual polar cap, dynamic dune fields, and so forth; (b) We make use of the topographic information from a reference DEM gridded from updated and adjusted MOLA profiles (Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022), and examine all available summer images that overlap with indicated sloped surfaces (with slope greater than an empirical threshold of 15°). If there are
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Locations of scarps with presence of ice blocks as examined in the High Resolution Imaging Science Experiment (HiRISE) images. Locations with slope greater than 15°, represented by dark orange points, are draped over the 1 km/pixel reference Digital Elevation Model (DEM) from reprocessed and then self-registered Mars Orbiter Laser Altimeter (MOLA) profiles (Xiao, Stark, Schmidt, et al., 2021). Images with the presence of ice blocks at locations where at least five repeated observations have been made by the HiRISE camera are shown with black background. Meanwhile, images examined to feature ice blocks along troughs and scarps indicated by MOLA slope map are exhibited with pink background. Images enclosed by white rectangles are taken as examples and enlarged in the corresponding insets to show the present ice blocks (Insets 1 and 2) and rocks (Insets 3 and 4), respectively. The scale shown in Inset 1 also applies to the other three insets. Images enclosed by red rectangles are examples with the presence of rocks in the circumpolar regions, including the polar erg, that is, massive dune fields surrounding the North Polar Layered Deposits (NPLD). The barchanoid dune field with rocks at (75°N, 282°E), studied in Mount and Titus (2015) to infer the depth of the seasonal deposits, is marked. Names of regions that feature clusterings of scarps with ice blocks are annotated. The Udzha Crater is marked of which only the topmost sharp-edged rims rise above the polar layered deposits to hint at its circular shape. Scarp 1 studied in this paper is also marked for reference.</figcaption></figure>

visible ice blocks, then we mark the corresponding scarp as positive. We set no height constraints on the ice blocks as they should be at least tall enough to resolve the thin layer of seasonal deposits during middle-to-late northern spring. Actually, the pixel size of HiRISE (down to 0.25 m) limits the detection of ice blocks and shadows under 0.71 m in size and length, respectively (2√2 ≈ 2.83 pixels). Thus, the images themselves have already filtered out small ice blocks. In reality, there should exist a large quantity of these small blocks, as indicated by the frequency-size distribution of detected newly fallen ice blocks, which generally follows a power law (Fanara et al., 2020b).

The reference DEM is derived based on the co-registration technique (Gläser et al., 2013; Stark et al., 2015). We apply the concept of “self-registration” to improve the positioning of individual laser profiles. A random subset of laser profiles (fixed at 0.25) is selected and then co-registered to a footprint point cloud formed from the remaining profiles. After enough repeats (set to 30), we effectively remove all offsets between the profiles due to residual errors in laser alignment calibration, spacecraft attitude, timing, and trajectory. In Text S3 in Supporting Information S1, we perform a probabilistic analysis on the number of iterations needed to be performed in the self-registration process which justifies the repeats for up to 30 times. The DEM is gridded with 1 km in pixel size and cannot fully resolve steep scarps (normally with a width of some hundreds of meters in the images or less depending on the steepness), instead the indicated slopes are average representative over large-scale topography, for example, the transition zone between the NPLD and the Basal Unit. The DEM pixels with slope greater than an empirical threshold of 15° are shown in Figure 15. This relatively low threshold is set to ensure that all scarps with
ice block presence will lie in the search path. The extracted sloped surfaces include the marginal scarps, with heights that range from 200 m to more than 1,200 m, delineated by Massé et al. (2012). Indeed, this relatively low threshold also leads to plenty of instances of spiral troughs being indicated as candidates to search for ice blocks.

We locate 138 sites of which there are, at least, five different observations and an additional of 210 sites, a majority of which feature no repeated observation or with only one extra coverage, along the sloped surfaces indicated by the MOLA terrain model. The summer images at those sites are manually examined and the ones with presence of ice blocks are shown in Figure 15. A total of 66 images over various locations show signs of recent activity (39 out of 138 and 27 out of 210 from the two categories, respectively). Active scarps with more than or equal to 5 repeated observations cluster in Olympia Rupes (including Scarp 1), Abalos Scopuli, and Boreum Cavus together with Tenius Cavus within Chasma Boreale. Apart from Olympia Rupes and Chasma Boreale, active scarps located through slope indicator extend to Rupes Tenuis and Gemini Scopuli. However, it is worth mentioning that the scarps at Olympia Rupes, Abalos Scopuli, and Chasma Boreale are much more dynamic than those at the other places, as indicated by the large number of aprons of debris and ice blocks spotted there. These mass-wasting activities are spatially correlated with Basal Unit exposure at their base (compare to Figure 2 in Nerozzi et al. (2022)). In contrast, scarps along Rupes Tenuis lack apparent aprons of debris but ice blocks reside right at the layered terrain of the scarps (refer to Inset 1 in Figure 15 for an example). Actually, Rupes Tenuis itself belongs to the exposed Basal Unit, and not the NPLD. However, to prevent unnecessary confusion, we do not attempt to distinguish these scarps from those of the NPLD. Scarps in Gemini Scopuli over Gemina Lingula also lack aprons of debris, and ice blocks are large in quantity and extremely uniformly distributed at close vicinity of the scarps (refer to Inset 2 in Figure 15 for an example). Comparing the summer images in MY29 and MY36 at the location of Inset 2 shows no apparent fracture-caused detachment of fragments over the scarp, indicating these ice blocks, being contiguous to the scarp, may be anteriorly emplaced. The scarps with visible ice blocks spatially cover latitudes from 78°N to 86°N. In terms of longitude, there exist a major gap from 27°E to 117°E over the Gemini Scopuli where Udzha Crater (81.8°N, 75.0°E) is the only site that has been spotted with ice blocks (refer to Figure 15 for location). It should be noted that the complete set of HiRISE images that have been examined do not cover all of the scarps. When searching through the images, only the ones acquired in the summertime have been downloaded and inspected. Additionally, there exist gaps in HiRISE coverage of the scarps which can be expected to be gradually bridged as more upcoming observations have been planned. As such, the set of scarps with presence of ice blocks presented in Figure 15 should stand as the lower limit on the maximum spatial extent of the expected outputs by the proposed approaches.

With ice blocks being spatially limited to high latitudes where the thickest seasonal layer is expected (Gary-Bicas et al., 2020; Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022), rocks at regions surrounding the north polar cap, either over the crater rims or between the dune fields in the polar erg, and at lower latitudes can also be utilized for the purpose of measuring the thickness of the seasonal layer (Cull et al., 2010; Mount & Titus, 2015). Some typical rock sites are enclosed by red rectangles in Figure 15, with Insets 3 and 4 showing zoom-in views of the rocks over the crater rim and between the dune fields, respectively. For large boulders over the craters and their vicinity, they could be ejected during the impacts. The dune field, named Dunes, within which the rocks have been utilized in Mount and Titus (2015) to infer the thickness of the seasonal deposits there is also taken as an example. Unfortunately, no rocks have been spotted between the massive linear dune fields in Olympia Undae where maximum depth up to ~4 m and off-season thickness variations up to ~3 m in magnitude have been claimed after examining the dynamic MOLA height records (Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022). Instead, bright materials between the dunes which can persist throughout the northern summer have been spotted. Existing literature suggested that they are most likely to be composed of gypsum and other evaporitic minerals as well (Fishbaugh et al., 2007; Horgan et al., 2009).

By applying the THB approach to all available sites with ice blocks and rocks in the polar region, we can obtain good spatial samplings of the thickness evolution of the seasonal deposits and shed comprehensive insights into the long-term vertical growth and retreat of the SNPC. In particular, these results can be used to examine whether the seasonal deposits are much shallower than the MOLA-derived results during springtime over the entire Martian north polar region (Section 5.5). If so, that means the average bulk density of the seasonal deposits in spring should be much higher than thought, which can have implications for snow/ice metamorphism and translucency of the deposits (e.g., Eluszkiewicz et al., 2005; Matsuo & Heki, 2009; Mount & Titus, 2015; Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022). The expected measurements can also shed light on the spatially variable proportions of direct condensation as frosts and atmospheric precipitation as snowfalls in forming the
seasonal layer. These proportions can help determine if CO₂ snowing storms are much more frequent and violent than thought (Section 5.6). Meanwhile, THV can be applied to specific (relatively tall) ice blocks to look into the thickness evolution at individual spots of interest. Possible significant interannual depth variations with magnitudes of more than 0.2 m can be identified and interpreted in the context of Martian climate change. Those enhanced understandings can place crucial constraints on the Martian volatile cycles and climate models (e.g., Forget et al., 2013; Mischna et al., 2003). Meanwhile, the thickness of the seasonal layer can assist in designing of future landers, rovers, or helicopters that are to drill in the spiral troughs and decipher the Late Amazonian climate of Mars stored in the NPLD (I. B. Smith et al., 2020; Matthies et al., 2022). As such, we strongly advocate for more future observations of the active scarps and circumpolar rock fields at the Martian north polar region using the HiRISE camera, including stereo pairs.

## 5.8. Limitations and Outlook

We would like to stress that THB can be applied to all HiRISE images acquired from MY29 on and at all of the active scarps across the north polar region to form thickness evolution time series in multiple years. Unfortunately, that would require a tremendous amount of efforts as it relies on manual identification of the proper bounding ice blocks and measurements of their shadows in the summer images. As such, we plan to build on previously developed automatic software for this purpose, by combining computer vision for delineating the shape of the ice blocks (Fanara et al., 2020a, 2020b) with Convolutional Neural Networks for identifying newly placed ice blocks (Su, Fanara, Xiao, et al., 2023). In order to infer their corresponding heights, extracted horizontal boundary of the ice blocks can be fitted using ellipsoids, triangular prisms, circular cones, or pyramids according to the delineated shape of the shadows. These software programs are established and well maintained in Institute of Geodesy and Geoinformation Science, Technische Universität Berlin. Existing tools that automatically locate and measure boulders on Mars (e.g., Golombek et al., 2008; Hood et al., 2022; Nagle-McNaughton et al., 2020), will be adapted for detecting ice blocks and used for cross-validation purposes. The main idea is to automatically identify the ice blocks, exclude newly formed ice blocks, and determine the corrected ratio of the ice blocks that have not been completely submerged during late winter and spring as compared to that in the summer. Meanwhile, the approximated sizes of the ice blocks, including approximated vertical dimension, during the summertime should be automatically determined. Then, we can statistically relate the detected ratio to the thickness of the seasonal ice layer according to probability distribution of the ice block heights:

$$ \int_{\Omega_Y}^{H_{\max}} p(H) dH = \Delta, \quad \text{s.t.} \quad \int_{0}^{H_{\max}} p(H) dH = 1, \qquad (20) $$

where $\Omega_Y$ is the thickness of the seasonal layer at $L_s = \Upsilon^\circ$ corresponding to the acquisition time of a specific HiRISE image, $H_{\max}$ is the maximum height of the ice blocks within an examined region, and $p(H)$ is the normalized probability density of the ice block heights. Meanwhile, $\Delta$ is the software-determined quantity ratio of ice blocks that have not been fully submerged at $L_s = \Upsilon^\circ$. It should be noted that Equation 20 should be in its discretized form as the number of ice blocks within an interest region does not reach infinity. The expected advantage of this automatic approach, tentatively termed automatic THB, is the capability to be efficiently applied to entire set of active scarps across the north polar region (refer to Section 5.7 for locations) and to monitor the depth evolution in real time once a HiRISE image beams back to Earth. An additional merit would be that it can completely avoid the cumbersome procedure of measuring dimensions of the shadows, as required in THV and THB.

Spatially and temporally inhomogeneous growth and retreat of the seasonal polar caps have been identified and analyzed over various seasons using the MOLA records (Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022; Xiao, Stark, Schmidt, Hao, Su, et al., 2022). Unfortunately, a realistic error propagation is not possible without ground truth for assessment. As a result, the unexpected features observed, for example, exceptionally high maximum depth (~4 m) and off-season depth trends observed within Olympia Undae, remain to be verified from independent dataset. Indeed, current accuracy assessment of the MOLA results rely merely on indirect evidence, for example, cross-validation of results from various teams or deviations from expected values (normally zero) at specific regions and solar longitudes (D. E. Smith et al., 2001; Aharonson et al., 2004; Xiao, Stark, Schmidt, Hao, Su, et al., 2022; Xiao, Stark, Schmidt, Hao, Steinbrügge, et al., 2022). The large biases in MOLA-derived results during spring as stated in Section 5.5 stress the importance of proper calibration of the MOLA results.
Although temporally limited to late winter and spring (Section 4) and spatially confined to high latitudes (Sec-
tion 5.7), the results by inspecting the ice blocks and rocks over the entire north polar region in HiRISE images
can serve to calibrate the existing MOLA results. Then, extrapolation of the calibrated MOLA results to un-
covered latitudes can be implemented in combination of boundary constraints, for example, CROCUS dates from
optical cameras and thermal infrared spectrometers (e.g., Schmidt et al., 2009; Piqueux et al., 2015). The term
introduced by Kieffer et al. (2000), is a helpful mnemonic marker, which identifies the time in spring at which the
snow/ice has completely disappeared. Meanwhile, these expected results can be utilized to validate the processing
procedure of the SHARAD radar sounding and altimetric measurements for the purpose of looking into the depth
evolution of the seasonal deposits. These radar records harbor the potential to decipher the long-term thickness
evolution patterns of the seasonal deposits at both poles of Mars (e.g., Raguso & Nunes, 2021; Steinbrügge
et al., 2021).

The depth evolution curves from THV and THB can achieve an uncertainty of better than 0.1 m at best scenarios, which is mainly limited by the spatial resolution of the HiRISE camera (down to 0.25 m). This dictates that delicate interannual variations, for example, less than ~0.2 m, cannot be confidently detected at the current stage. Thus, future super-resolution cameras to orbit Mars would dramatically enhance the performance of the proposed approaches to delve into the snow/ice depth evolution patterns. Indeed, the Mars Next Orbiter Science Analysis Group (NEX-SAG) already identified a baseline instrument of a super-high-resolution optical imager (HiRISE class at 30 cm/pixel; or even with a better resolution of 10–15 cm/pixel) to reveal detailed morphology over limited areas for science and site reconnaissance (Zurek et al., 2021). This recommended imager is capable of, for example, capturing stratigraphic details, and thus the Late Amazonian climate records, stored in the PLD (I. B. Smith et al., 2020) and better observing various surface dynamics such as those related to slopes (Dundas et al., 2021).

Using MOLA-measured thickness variations of the seasonal deposits from Xiao, Stark, Schmidt, Hao, Su, et al. (2022) and Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2022), Wagner et al. (2023) modeled the loading-induced lithospheric deflection based on mechanical properties of the Martian interior. They found that the deflection can be up to ~6 cm in magnitude at the poles and can be unambiguously detected by upcoming Interferometric Synthetic Aperture Radars (InSARs) or next-generation laser altimeters. We note that the thickness variations of the seasonal deposits from MOLA altimetric records were retrieved by examining the surface height changes, which actually represent the combined effects of CO₂ deposition/sublimation and the resultant lithospheric deflection. In contrast, thickness measurements from THV and THB are obtained by examining the shadows of the ice blocks and thus remain unaffected by the lithospheric deflection. It is anticipated that we can extract the periodic lithospheric deflection by comparing the surface height variations from future InSARs or laser altimeters with the seasonal deposit thicknesses from examining the ice blocks. Constraints on the lithospheric deflection can inform us about the structure and rheology of the Martian crust and mantle which have important implications for Mars' thermal history. It will be challenging to differentiate between currently plausible Mars interior models as the precision in measuring the deflection would need to be in the millimeter or even sub-millimeter level (Wagner et al., 2023). However, this proposed concept does harbor some potential by including a large amount of independent surface height variation and snow/ice depth measurements for spatial and temporal averaging purposes. The fact that the lithospheric deflection peaks at high polar latitudes where ice blocks are located further enhances the feasibility of the measuring concept. Meanwhile, a future super-high-resolution optical imager to Mars would greatly build up the prospects (Zurek et al., 2021). A factor that will need to be considered and thoroughly modeled in the determination of the subtle lithospheric deflection is the solid body tidal deformation from the Sun and Phobos which maximizes at around the equator and reduces to less than 1 cm at the poles (Wagner et al., 2023).

**6. Conclusion**

We use the shadow variations of the ice blocks at the foot of the NPLD steep scarps, complementing that of the rocks, to infer the vertical evolution of the Martian seasonal deposits. We relate the shadow length of an ice block to its height using a rigorous geometric model, which is based on bundle-adjusted and orthorectified high-resolution HiRISE images and takes both the solar and slope properties into consideration. Building on this model, we present two independent and yet complementary approaches: (a) THV that subtracts the ice block heights in the summer to that observed in the spring; (b) THB that locates ice blocks that have been completely covered to place lower limits on the thickness of the seasonal deposits, and ice blocks that still stick out of the
seasonal cover to put upper limits. Our interpretation is based on reasonable assumptions about the distribution of snowfall and frost around the ice blocks and their surroundings. However, it should be cautioned that deviations from these ideal circumstances are possible and future work is needed to quantify their impacts. We experimentally apply the methods to a steep scarp centered at (85.0°N, 151.5°E). The results show that the average thickness due to accumulation of snowfalls in MY31 is 0.97 ± 0.13 m at Ls = 350.7° in late winter, 0.64 ± 0.08 m at Ls = 7.0° in early spring, 0.21 ± 0.05 m at Ls = 42.8° in middle spring, and gradually decreases toward summer solstice (Ls = 90°). The maximum snowfall thickness of ~1 m is two orders of magnitude larger than theoretical values predicted by Alsaeed and Hayne (2022). We proceed to use the widening of the ice blocks as a proxy to approximate the depth of the directly condensed layer over the ice blocks. We show that average thickness of the frost layer in MY31 reaches 0.64 ± 0.18 m at Ls = 350.7° in late winter, quasi-linearly declines to 0.26 ± 0.03 m at Ls = 42.8° in middle spring and to 0.045 ± 0.035 m at Ls = 69.5° in late spring. We thus show that atmospheric deposition as snowfalls, compared to direct surface condensation as frosts, can even contribute more to the thickness and volume of the seasonal snow/ice layer at Scarp 1 during wintertime. The aggregate thickness of the seasonal deposits at Scarp 1 in MY31 then stands at 1.63 ± 0.22 m at Ls = 350.7° in late winter which then gradually declines to 0.45 ± 0.06 m at Ls = 42.8° in middle spring and 0.06 ± 0.05 m at Ls = 69.5° in late spring. The majority of the associated uncertainties can be better than ~0.1 m. Compared to the existing MOLA results, our thickness estimates are consistently lower, by a magnitude of up to 0.8 m, throughout the spring. We attribute a large portion of these inconsistencies to biases inherited in the MOLA-derived thickness measurements. The available HiRISE images can temporally span from MY29 to MY36 (2008–2021), and we observe the snowfall thickness in the very early spring in MY36 is 1.01 ± 0.10 m, exceeding that in MY31 by a magnitude of 0.36 ± 0.13 m.

Note that the results in this study are obtained at one single scarp at the Martian north polar region. Building on this study, we will develop automatic version of the approaches and extend the measurements to all scarps with presence of ice blocks, covered by multiple HiRISE images at the same time, at the Martian north polar region. We will also apply the proposed approaches to the circumpolar regions with presence of rocks to extend the thickness evolution measurements to lower polar latitudes. These expected results can answer the question that if the seasonal cap over the entire Martian north polar region during its recession in spring is much shallower than we previously thought. Meanwhile, the proportions of snowfall and frosts in constituting the seasonal deposits over the entire polar region can be expected to be revealed. These future results can also serve to calibrate previous MOLA results and validate contemporary anticipated SHARAD results. Besides, these expected outcomes can help to examine the long-term spatial-temporal heterogeneity of the surface-atmosphere exchange and determine the mass balance of the NPLD. As such, we suggest more frequent targeting of the HiRISE camera to the active scarps where ice blocks exist and circumpolar fields with presence of rocks. Meanwhile, anticipated super-high-resolution optical cameras to Mars would enable unambiguous identification of small interannual depth variations, for example, less than 0.2 m in magnitude, of the seasonal deposits. Thickness measurements from examining the ice blocks, in combination with high-precision surface height variations from future InSARs or laser altimeters, have the potential of inferring the polar dynamic lithospheric deflection which can shed light on the Martian interior.

## Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

## Data Availability Statement

HiRISE images are available at NASA Planetary Data System (PDS) (McEwen, 2007). The generated HiRISE DEM and scripts to reproduce Figures 9, 11, 13 and 14 are published in Mendeley Data (Xiao et al., 2024). The self-registered MOLA reference DEM and thickness variation time series from Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2022) at the Martian north polar region are published in Mendeley Data (Xiao, Stark, Schmidt, et al., 2021). The SPICE kernels used for calculation of the solar condition are available at NASA PDS (Semenov et al., 1998). The MGS MOC global-map images can be obtained at the MOC Image Gallery maintained by Malin Space Science Systems (Malin et al., 2001). To access these images, navigate to https://asimov.msss.com/moc_ gallery/, and select the category of "MOC Wide Angle Images Regional Views and Context Images".
## Acknowledgments

The authors would like to thank Nicholas L. Wagner (Baylor University), Noora R. Alsaeed (University of Colorado Boulder), and Lida Fanara (DLR) for helpful discussions. We thank Anni Määttänen for handling the manuscript, Mike Sori (Purdue University) and two anonymous reviewers for constructive comments and suggestions which helped to significantly improve the manuscript. We also acknowledge the MRO HiRISE and MGS MOC teams for providing the great products that enabled this research. HX, LML, and PJG acknowledge financial support from project PID2021-126365NB-C21 (MCI/AEI/FEDER, UE) and from the grant CEX2021-001131-S funded by MCI/AEI/10.13039/501100011033. SS thanks the China Scholarship Council (CSC) for financial support of studying in Germany. FS acknowledges support from the "Institut National des Sciences de l'Univers" (INSU), the "Centre National de la Recherche Scientifique" (CNRS), and "Centre National d'Etudes Spatiales" (CNES) through the "Programme National de Planétologie." Open Access funding enabled and organized by Projekt DEAL.

## References

Acton, C. H. J. (1996). Ancillary data services of NASA's navigation and ancillary information facility. *Planetary and Space Science*, 44(1), 65–70. https://doi.org/10.1016/0032-0633(95)00107-7

Aharonson, O., Zuber, M. T., Smith, D. E., Neumann, G. A., Feldman, W. C., & Prettyman, T. H. (2004). Depth, distribution, and density of CO₂ deposition on Mars. *Journal of Geophysical Research: Planets*, 109(E5). https://doi.org/10.1029/2003JE002223

Alsaeed, N. R., & Hayne, P. O. (2022). Transport of water into the polar regions of Mars through scavenging by CO₂ snowfall. *Journal of Geophysical Research: Planets*, 127(11), e2022JE007386. https://doi.org/10.1029/2022JE007386

Andrieu, F., Schmidt, F., Douté, S., & Chassefière, E. (2018). Ice state evolution during spring in Richardson crater, Mars. *Icarus*, 315, 158–173. https://doi.org/10.1016/j.icarus.2018.06.019

Appéré, T., Schmitt, B., Langevin, Y., Douté, S., Pommerol, A., Forget, F., et al. (2011). Winter and spring evolution of northern seasonal deposits on Mars from OMEGA on Mars Express. *Journal of Geophysical Research: Planets*, 116(E5), E50001. https://doi.org/10.1029/2010JE003762

Becerra, P., Smith, I. B., Hibbard, S., Andres, C., Bapst, J., Bramson, A. M., et al. (2021). Past, present, and future of Mars polar science: Outcomes and outlook from the 7th International Conference on Mars Polar Science and Exploration. *Planetary Science Journal*, 2(5), 209. https://doi.org/10.3847/PSJ/ac19a5

Becerra, P., Sori, M. M., & Byrne, S. (2017). Signals of astronomical climate forcing in the exposure topography of the North Polar Layered Deposits of Mars. *Geophysical Research Letters*, 44(1), 62–70. https://doi.org/10.1002/2016GL071197

Benson, J. L., Kass, D. M., & Kleinböhl, A. (2011). Mars' north polar hood as observed by the Mars Climate Sounder. *Journal of Geophysical Research: Planets*, 116(E3), E03008. https://doi.org/10.1029/2010JE003693

Beyer, R. A., Alexandrov, O., & McMichael, S. (2018). The Ames Stereo Pipeline: NASA's open source software for deriving and processing terrain data. *Earth and Space Science*, 5(9), 537–548. https://doi.org/10.1029/2018EA000409

Blackburn, D. G., Bryson, K. L., Chevrier, V. F., Roe, L. A., & White, K. F. (2010). Sublimation kinetics of CO₂ ice on Mars. *Planetary and Space Science*, 58(5), 780–791. https://doi.org/10.1016/j.pss.2009.12.004

Botev, Z. I., Grotowski, J. F., & Kroese, D. P. (2010). Kernel density estimation via diffusion. *Annals of Statistics*, 38(5), 2916–2957. https://doi.org/10.1214/10-AOS799

Brothers, S. C., Kocurek, G., & Holt, J. W. (2018). Sequence architecture of the cavi unit, Chasma Boreale, Mars. *Icarus*, 308, 42–60. https://doi.org/10.1016/j.icarus.2017.06.024

Brothers, T. C., Holt, J. W., & Spiga, A. (2015). Planum Boreum basal unit topography, Mars: Irregularities and insights from SHARAD. *Journal of Geophysical Research: Planets*, 120(7), 1357–1375. https://doi.org/10.1002/2015JE004830

Brown, A. J., Calvin, W. M., Becerra, P., & Byrne, S. (2016). Martian north polar cap summer water cycle. *Icarus*, 277, 401–415. https://doi.org/10.1016/j.icarus.2016.05.007

Byrne, S., Sori, M. M., Russell, P. S., Pathare, A. V., Becerra, P., Molaro, J. L., et al. (2017). Mars polar cliffs: Stressed out and falling apart. In *European planetary science congress* (Vol. 11).

Calvin, W. M., James, P. B., Cantor, B. A., & Dixon, E. M. (2015). Interannual and seasonal changes in the north polar ice deposits of Mars: Observations from MY 29–31 using MARCI. *Icarus*, 251, 181–190. https://doi.org/10.1016/j.icarus.2014.08.026

Colaprete, A., Barnes, J. R., Haberle, R. M., Hollingsworth, J. L., Kieffer, H. H., & Titus, T. N. (2005). Albedo of the south pole on Mars determined by topographic forcing of atmosphere dynamics. *Nature*, 435(7039), 184–188. https://doi.org/10.1038/nature03561

Colaprete, A., Barnes, J. R., Haberle, R. M., & Montmessin, F. (2008). CO₂ clouds, CAPE and convection on Mars: Observations and general circulation modeling. *Planetary and Space Science*, 56(2), 150–180. https://doi.org/10.1016/j.pss.2007.08.010

Cull, S., Arvidson, R. E., Mellon, M., Wiseman, S., Clark, R., Titus, T., et al. (2010). Seasonal H₂O and CO₂ ice cycles at the Mars Phoenix landing site: 1. Prelanding CRISM and HiRISE observations. *Journal of Geophysical Research: Planets*, 115(E4). https://doi.org/10.1029/2009JE003340

Defraigne, P., de Viron, O., Dehant, V., Van Hoolst, T., & Hourdin, F. (2000). Mars rotation variations induced by atmosphere and ice caps. *Journal of Geophysical Research: Planets*, 105(E10), 24563–24570. https://doi.org/10.1029/1999JE001227

Dundas, C. M., Becerra, P., Byrne, S., Chojnacki, M., Daubar, I. J., Diniega, S., et al. (2021). Active Mars: A dynamic world. *Journal of Geophysical Research: Planets*, 126(8), e2021JE006876. https://doi.org/10.1029/2021JE006876

Eluszkiewicz, J., Moncet, J.-L., Titus, T. N., & Hansen, G. B. (2005). A microphysically-based approach to modeling emissivity and albedo of the Martian seasonal caps. *Icarus*, 174(2), 524–534. https://doi.org/10.1016/j.icarus.2004.05.025

Fanara, L., Gwinner, K., Hauber, E., & Oberst, J. (2020a). Automated detection of block falls in the north polar region of Mars. *Planetary and Space Science*, 180, 104733. https://doi.org/10.1016/j.pss.2019.10.4733

Fanara, L., Gwinner, K., Hauber, E., & Oberst, J. (2020b). Present-day erosion rate of north polar scarps on Mars due to active mass wasting. *Icarus*, 342, 113434. https://doi.org/10.1016/j.icarus.2019.113434

Fishbaugh, K. E., & Head III, J. W. (2005). Origin and characteristics of the Mars north polar basal unit and implications for polar geologic history. *Icarus*, 174(2), 444–474. https://doi.org/10.1016/j.icarus.2004.06.021

Fishbaugh, K. E., Poulet, F., Chevrier, V., Langevin, Y., & Bibring, J.-P. (2007). On the origin of gypsum in the Mars north polar region. *Journal of Geophysical Research: Planets*, 112(E7). https://doi.org/10.1029/2006JE002862

Forget, F., Wordsworth, R., Millour, E., Madeleine, J.-B., Kerber, L., Leconte, J., et al. (2013). 3D modelling of the early Martian climate under a denser CO₂ atmosphere: Temperatures and CO₂ ice clouds. *Icarus*, 222(1), 81–99. https://doi.org/10.1016/j.icarus.2012.10.019

Gary-Bicas, C. E., Hayne, P. O., Horvath, T., Heavens, N. G., Kass, D. M., Kleinböhl, A., et al. (2020). Asymmetries in snowfall, emissivity, and albedo of Mars' seasonal polar caps: Mars climate sounder observations. *Journal of Geophysical Research: Planets*, 125(5), e2019JE006150. https://doi.org/10.1029/2019JE006150

Gläser, P., Haase, I., Oberst, J., & Neumann, G. A. (2013). Co-registration of laser altimeter tracks with digital terrain models and applications in planetary science. *Planetary and Space Science*, 89, 111–117. https://doi.org/10.1016/j.pss.2013.09.012

Golombek, M. P., Huertas, A., Marlow, J., Mcgrane, B., Klein, C., Martinez, M., et al. (2008). Size-frequency distributions of rocks on the northern plains of Mars with special reference to Phoenix landing surfaces. *Journal of Geophysical Research: Planets*, 113(E3). https://doi.org/10.1029/2007JE003065

Grima, C., Kofman, W., Mouginot, J., Phillips, R. J., Hérique, A., Biccari, D., et al. (2009). North polar deposits of Mars: Extreme purity of the water ice. *Geophysical Research Letters*, 36(3). https://doi.org/10.1029/2008GL036326

Grisolle, F. (2013). Seasonal condensates on Mars: Experimental study of formation and metamorphism of CO₂ ices. PhD thesis, Université de Grenoble. Retrieved from https://tel.archives-ouvertes.fr/tel-01010519

Hansen, C., Byrne, S., Portyankina, G., Bourke, M., Dundas, C., McEwen, A., et al. (2013). Observations of the northern seasonal polar cap on Mars: I. Spring sublimation activity and processes. *Icarus*, 225(2), 881–897. https://doi.org/10.1016/j.icarus.2012.09.024
Hansen, C., Thomas, N., Portyankina, G., McEwen, A., Becker, T., Byrne, S., et al. (2010). HiRISE observations of gas sublimation-driven activity in Mars' southern polar regions: I. Erosion of the surface. *Icarus*, *205*(1), 283–295. https://doi.org/10.1016/j.icarus.2009.07.021

Hayne, P. O., Paige, D. A., & Heavens, N. G., & the Mars Climate Sounder Science Team. (2014). The role of snowfall in forming the seasonal ice caps of Mars: Models and constraints from the Mars Climate Sounder. *Icarus*, *231*, 122–130. https://doi.org/10.1016/j.icarus.2013.10.020

Hayne, P. O., Paige, D. A., Schofield, J. T., Kass, D. M., Kleinböhl, A., Heavens, N. G., & McCleese, D. J. (2012). Carbon dioxide snow clouds on Mars: South polar winter observations by the Mars Climate Sounder. *Journal of Geophysical Research: Planets*, *117*(E8). https://doi.org/10.1029/2011JE004040

Hayward, R. K., Fenton, L. K., Tanaka, K. L., Titus, T. N., Colaprete, A., & Christensen, P. R. (2010). Mars global digital dune database: *MC-1*, *US Geological Survey open-file report 2010-1170*. US Geological Survey. https://doi.org/10.3133/ofr20101170

Head III, J. W., Hiesinger, H., Ivanov, M. A., Kreslavsky, M. A., Pratt, S., & Thomson, B. J. (1999). Possible ancient oceans on Mars: Evidence from Mars Orbiter Laser Altimeter data. *Science*, *286*(5447), 2134–2137. https://doi.org/10.1126/science.286.5447.2134

Hepburn, A. J., Holt, T., Hubbard, B., & Ng, F. (2019). Creating HiRISE digital elevation models for Mars using the open-source Ames Stereo Pipeline. *Geoscientific Instrumentation, Methods and Data Systems*, *8*(2), 293–313. https://doi.org/10.5194/gi-8-293-2019

Hood, D. R., Sholes, S. F., Karunatillake, S., Fassett, C. I., Ewing, R. C., & Levy, J. (2022). The Martian boulder automatic recognition system, MBARS. *Earth and Space Science*, *9*(9), e2022EA002410. https://doi.org/10.1029/2022EA002410

Horgan, B. H., Bell III, J. F., Noe Dobrea, E. Z., Cloutis, E. A., Bailey, D. T., Craig, M. A., et al. (2009). Distribution of hydrated minerals in the north polar region of Mars. *Journal of Geophysical Research: Planets*, *114*(E1). https://doi.org/10.1029/2008JE003187

Kieffer, H. H., Titus, T. N., Mullins, K. F., & Christensen, P. R. (2000). Mars south polar spring and summer behavior observed by TES: Seasonal cap evolution controlled by frost grain size. *Journal of Geophysical Research: Planets*, *105*(E4), 9653–9699. https://doi.org/10.1029/1999JE001136

Laskar, J., Levrard, B., & Mustard, J. F. (2002). Orbital forcing of the Martian Polar Layered Deposits. *Nature*, *419*(6905), 375–377. https://doi.org/10.1038/nature01066

Leighton, R. B., & Murray, B. C. (1966). Behavior of carbon dioxide and other volatiles on Mars. *Science*, *153*(3732), 136–144. https://doi.org/10.1126/science.153.3732.136

Le Maistre, S., Rivoldini, A., Caldiero, A., Yseboodt, M., Baland, R.-M., Beuthe, M., et al. (2023). Spin state and deep interior structure of Mars from InSight radio tracking. *Nature*, *619*(7971), 1–5. https://doi.org/10.1038/s41586-023-06150-0

Leys, C., Ley, C., Klein, O., Bernard, P., & Licata, L. (2013). Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median. *Journal of Experimental Social Psychology*, *49*(4), 764–766. https://doi.org/10.1016/j.jesp.2013.03.013

Määttänen, A., & Montmessin, F. (2021). Clouds in the Martian atmosphere. In *Oxford research encyclopedia of planetary science*. https://doi.org/10.1093/acrefore/9780190647926.013.114

Malin, M. C., Bell, J. F., Cantor, B. A., Caplinger, M. A., Calvin, W. M., Clancy, R. T., et al. (2007). Context Camera investigation on board the Mars Reconnaissance Orbiter. *Journal of Geophysical Research: Planets*, *112*(E5). https://doi.org/10.1029/2006JE002808

Malin, M. C., Edgett, K. S., Davis, S. D., Caplinger, M. A., Jensen, E., Supulver, K. D., et al. (2001). Malin Space Science Systems Mars Orbiter Camera image gallery [Dataset]. *Malin Space Science Systems*. http://www.mss.com/moc gallery/2001

Massé, M., Bourgeois, O., Le Mouélic, S., Verpoorter, C., Spiga, A., & Le Deit, L. (2012). Wide distribution and glacial origin of polar gypsum on Mars. *Earth and Planetary Science Letters*, *317*, 44–55. https://doi.org/10.1016/j.epsl.2011.11.035

Matsuo, K., & Heki, K. (2009). Seasonal and inter-annual changes of volume density of Martian CO₂ snow from time-variable elevation and gravity. *Icarus*, *202*(1), 90–94. https://doi.org/10.1016/j.icarus.2009.02.023

Matthies, L., Kennett, A., Kerber, L., Fraeman, A., & Anderson, R. C. (2022). Prospects for very long-range Mars rover missions. In *IEEE aerospace conference (AERO)* (pp. 1–11). https://doi.org/10.1109/AERO53065.2022.9843681

McEwen, A. S. (2007). Mars Reconnaissance Orbiter High Resolution Imaging Science Experiment, reduced data record, MRO-M-HIRISE-3-RDR-V1.1 [Dataset]. *NASA Planetary Data System*. https://doi.org/10.17189/1520303

McEwen, A. S., Eliason, E. M., Bergstrom, J. W., Bridges, N. T., Hansen, C. J., Delamere, W. A., et al. (2007). Mars Reconnaissance Orbiter's high resolution imaging science experiment (HiRISE). *Journal of Geophysical Research: Planets*, *112*(E5). https://doi.org/10.1029/2005JE002605

McGlasson, R. A., Bramson, A. M., Morgan, G. A., & Sori, M. M. (2023). Varied histories of outlier polar ice deposits on Mars. *Journal of Geophysical Research: Planets*, *128*(5), e2022JE007592. https://doi.org/10.1029/2022JE007592

Mc Keown, L. E., Diniega, S., Portyankina, G., Hansen, C. J., Aye, K. M., Piqueux, S., & Scully, J. E. C. (2023). Martian araneiforms: A review. *Journal of Geophysical Research: Planets*, *128*(4). https://doi.org/10.1029/2022JE007684

Métivier, L., Karatekin, Ö., & Dehant, V. (2008). The effect of the internal structure of Mars on its seasonal loading deformations. *Icarus*, *194*(2), 476–486. https://doi.org/10.1016/j.icarus.2007.12.001

Milkovich, S. M., Byrne, S., & Russell, P. S. (2012). Variations in surface texture of the north polar residual cap of Mars. In *Lunar and planetary science conference* (p. 2226).

Milkovich, S. M., & Head III, J. W. (2006). Surface textures of Mars' north polar layered deposits: A framework for interpretation and future exploration. *Mars*, *2*, 21–45. https://doi.org/10.1555/mars.2006.0003

Mischna, M. A., Richardson, M. I., Wilson, R. J., & McCleese, D. J. (2003). On the orbital forcing of Martian water and CO₂ cycles: A general circulation model study with simplified volatile schemes. *Journal of Geophysical Research: Planets*, *108*(E6). https://doi.org/10.1029/2003JE002051

Mount, C. P., & Titus, T. N. (2015). Evolution of Mars' northern polar seasonal CO₂ deposits: Variations in surface brightness and bulk density. *Journal of Geophysical Research: Planets*, *120*(7), 1252–1266. https://doi.org/10.1002/2014JE004706

Nagle-McNaughton, T. P., Williams, J. M., Gallegos, Z. E., Wilkie, H. A., Martinez, D. C., & Scuderi, L. A. (2020). Geographic information system based detection and quantification of boulders using HiRISE imagery: A case study in Jezero Crater. *Journal of Applied Remote Sensing*, *14*(1), 014522. https://doi.org/10.1117/1.JRS.14.014522

Navarro, T., Madeleine, J. B., Forget, F., Spiga, A., Millour, E., Montmessin, F., & Määttänen, A. (2014). Global climate modeling of the Martian water cycle with improved microphysics and radiatively active water ice clouds. *Journal of Geophysical Research: Planets*, *119*(7), 1479–1495. https://doi.org/10.1002/2013JE004550

Nerozzi, S., & Holt, J. W. (2019). Buried ice and sand caps at the north pole of Mars: Revealing a record of climate change in the cavi unit with SHARAD. *Geophysical Research Letters*, *46*(13), 7278–7286. https://doi.org/10.1029/2019GL082114

Nerozzi, S., Ortiz, M. R., & Holt, J. W. (2022). The north polar basal unit of Mars: An Amazonian record of surface processes and climate events. *Icarus*, *373*, 114716. https://doi.org/10.1016/j.icarus.2021.114716

Neukum, G., & Jaumann, R. (2004). HRSC: The high resolution stereo camera of Mars Express. *Mars Express: The Scientific Payload*, *1240*, 17–35.
Neumann, G. A., Abshire, J. B., Aharonson, O., Garvin, J. B., Sun, X., & Zuber, M. T. (2003). Mars Orbiter Laser Altimeter pulse width measurements and footprint-scale roughness. *Geophysical Research Letters*, 30*(11), 1561. https://doi.org/10.1029/2003GL017048

Neumann, G. A., Rowlands, D. D., Lemoine, F. G., Smith, D. E., & Zuber, M. T. (2001). Crossover analysis of Mars Orbiter Laser Altimeter data. *Journal of Geophysical Research: Planets*, 106*(E10), 23753–23768. https://doi.org/10.1029/2000JE001381

Ojha, L., Nerozzi, S., & Lewis, K. (2019). Compositional constraints on the north polar cap of Mars from gravity and topography. *Geophysical Research Letters*, 46*(15), 8671–8679. https://doi.org/10.1029/2019GL082294

Piqueux, S., Kleinböhl, A., Hayne, P. O., Kass, D. M., Schofield, J. T., & McCleese, D. J. (2015). Variability of the Martian seasonal CO₂ cap extent over eight Mars Years. *Icarus*, 251, 164–180. https://doi.org/10.1016/j.icarus.2014.10.045

Plaut, J. J., Picardi, G., Safaeinili, A., Ivanov, A. B., Milkovich, S. M., Cicchetti, A., et al. (2007). Subsurface radar sounding of the South Polar layered deposits of Mars. *Science*, 316*(5821), 92–95. https://doi.org/10.1126/science.1139672

Pommerol, A., Appéré, T., Portyankina, G., Aye, K.-M., Thomas, N., & Hansen, C. J. (2013). Observations of the northern seasonal polar cap on Mars III: CRISM/HiRISE observations of spring sublimation. *Icarus*, 225(2), 911–922. https://doi.org/10.1016/j.icarus.2012.08.039

Portyankina, G., Markiewicz, W. J., Thomas, N., Hansen, C. J., & Milazzo, M. (2010). HiRISE observations of gas sublimation-driven activity in Mars' southern polar regions: III. Models of processes involving translucent ice. *Icarus*, 205(1), 311–320. https://doi.org/10.1016/j.icarus.2009.08.029

Portyankina, G., Merrison, J., Iversen, J., Yoldi, Z., Hansen, C., Aye, K.-M., et al. (2019). Laboratory investigations of the physical state of CO₂ ice in a simulated Martian environment. *Icarus*, 322, 210–220. https://doi.org/10.1016/j.icarus.2018.04.021

Raguso, M. C., & Nunes, D. C. (2021). Monitoring of CO₂ seasonal variations at the Martian polar ice caps using SHARAD data. In *Lunar and planetary science conference* (p. 2438).

Russell, P. S., Byrne, S., & Hansen, C. J. (2010). Active mass wasting of ice layers and seasonal CO₂ frost in the north polar region of Mars. In *Lunar and planetary science conference*.

Russell, P. S., Byrne, S., Pathare, A., & Herkenhoff, K. E. (2012). Active erosion and evolution of Mars north polar scarps. In *Lunar and planetary science conference* (p. 2747).

Russell, P. S., Thomas, N., Byrne, S., Herkenhoff, K., Fishbaugh, K., Bridges, N., et al. (2008). Seasonally active frost-dust avalanches on a north polar scarp of Mars captured by HiRISE. *Geophysical Research Letters*, 35*(23). https://doi.org/10.1029/2008GL035790

Scheller, E., Ehlmann, B., Hu, R., Adams, D., & Yung, Y. (2021). Long-term drying of Mars by sequestration of ocean-scale volumes of water in the crust. *Science*, 372*(6537), 56–62. https://doi.org/10.1126/science.abc7717

Schmidt, F., Douté, S., Schmitt, B., Vincendon, M., Bibring, J.-P., & Langevin, Y., & the OMEGA Team. (2009). Albedo control of seasonal South Polar cap recession on Mars. *Icarus*, 200(2), 374–394. https://doi.org/10.1016/j.icarus.2008.12.014

Schmidt, F., & Portyankina, G. (2018). The exotic processes driving ephemeral seasonal surface change on Mars. In *Dynamic Mars: Recent and current landscape evolution of the red planet* (pp. 157–186). Elsevier. https://doi.org/10.1016/B978-0-12-813018-6.00005-4

Schmidt, F., Way, M. J., Costard, F., Bouley, S., Séjourné, A., & Aleinov, I. (2022). Circumpolar ocean stability on Mars 3 Gy ago. *Proceedings of the National Academy of Sciences*, 119*(4), e2112930118. https://doi.org/10.1073/pnas.2112930118

Semenov, B. V., Elson, L. S., & Acton, C. H. (1998). Mars global surveyor SPICE Kernels V1.0, MGS-M-SPICE-6-V1.0 [Dataset]. NASA *Planetary Data System*. https://doi.org/10.17189/1520101

Silverman, B. W. (1986). *Density estimation for statistics and data analysis* (Vol. 26). CRC Press.

Smith, D. E., Zuber, M. T., & Neumann, G. A. (2001). Seasonal variations of snow depth on Mars. *Science*, 294*(5549), 2141–2146. https://doi.org/10.1126/science.1066556

Smith, I. B., Hayne, P. O., Byrne, S., Becerra, P., Kahre, M., Calvin, W., et al. (2020). The Holy Grail: A road map for unlocking the climate record stored within Mars' polar layered deposits. *Planetary and Space Science*, 184(May 2019), 104841. https://doi.org/10.1016/j.pss.2020.104841

Smith, I. B., & Holt, J. W. (2010). Onset and migration of spiral troughs on Mars revealed by orbital radar. *Nature*, 465*(7297), 450–453. https://doi.org/10.1038/nature09049

Sori, M. M., Becerra, P., Bapst, J., Byrne, S., & McGlasson, R. A. (2022). Orbital forcing of Martian climate revealed in a south polar outlier ice deposit. *Geophysical Research Letters*, 49*(6), e2021GL097450. https://doi.org/10.1029/2021GL097450

Sori, M. M., Byrne, S., Hamilton, C. W., & Landis, M. E. (2015). Viscous flow rates of icy topography on the north polar layered deposits of Mars. *Geophysical Research Letters*, 43*(2), 541–549. https://doi.org/10.1002/2015GL067298

Stark, A., Oberst, J., Preusker, F., Gwinner, K., Peale, S. J., Margot, J.-L., et al. (2015). Mercury's rotational parameters from MESSENGER image and laser altimeter data: A feasibility study. *Planetary and Space Science*, 117, 64–72. https://doi.org/10.1016/j.pss.2015.05.006

Steinbrügge, G., Haynes, M. S., Schroeder, D. M., Scanlan, K. M., Stark, A., Young, D. A., et al. (2021). Altimetry measurements from planetary radar sounders and application to SHARAD on Mars. *IEEE Transactions on Geoscience and Remote Sensing*, 60, 1–14. https://doi.org/10.1109/TGRS.2021.3134638

Su, S., Fanara, L., Xiao, H., Hauber, E., & Oberst, J. (2023). Detection of detached ice-fragments at Martian polar scarps using a Convolutional Neural Network. *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, 16, 1728–1739. https://doi.org/10.1109/JSTARS.2023.3238968

Su, S., Fanara, L., Zhang, X., Hauber, E., Gwinner, K., & Oberst, J. (2023). Searching for the sources of ice block falls at the Martian north polar scarps. *Icarus*, 390, 115321. https://doi.org/10.1016/j.icarus.2022.115321

Thomas, N., Cremonese, G., Ziethe, R., Gerber, M., Brändli, M., Bruno, G., et al. (2017). The Colour and Stereo Surface Imaging System (CaSSIS) for the ExoMars Trace Gas Orbiter. *Space Science Reviews*, 212(3), 1897–1944. https://doi.org/10.1007/s11214-017-0421-1

Thomas, P. C., Calvin, W., Cantor, B., Haberle, R., James, P. B., & Lee, S. W. (2016). Mass balance of Mars' residual south polar cap from CTX images and other data. *Icarus*, 268, 118–130. https://doi.org/10.1016/j.icarus.2015.12.038

Van den Acker, E., Van Hoolst, T., de Viron, O., Defraigne, P., Forget, F., Hourdin, F., & Dehant, V. (2002). Influence of the seasonal winds and the CO₂ mass exchange between atmosphere and polar caps on Mars' rotation. *Journal of Geophysical Research: Planets*, 107*(E7), 9-1. https://doi.org/10.1029/2000JE001539

Wagner, N. L., James, P. B., Ermakov, A. E., & Sori, M. M. (2023). Evaluating the use of seasonal surface displacements and time-variable gravity to constrain the interior of Mars. *Authorea*. https://doi.org/10.22541/essoar.169230247.71360068/v1

Wang, H., & Ingersoll, A. P. (2002). Martian clouds observed by Mars Global Surveyor Mars Orbiter Camera. *Journal of Geophysical Research: Planets*, 107*(E10), 8-1. https://doi.org/10.1029/2001JE001815

Wilcoski, A. X., & Hayne, P. O. (2020). Surface roughness evolution and implications for the age of the north polar residual cap of Mars. *Journal of Geophysical Research: Planets*, 125(12), e2020JE006570. https://doi.org/10.1029/2020JE006570

Xiao, H., Stark, A., Schmidt, F., Hao, J., Steinbrügge, G., Cheng, Y., & Oberst, J. (2021). Spatio-temporal level variations of the Martian seasonal North Polar cap [Dataset]. *Mendeley Data*. https://doi.org/10.17632/x953mzxxvv.1
Xiao, H., Stark, A., Schmidt, F., Hao, J., Steinbrügge, G., Wagner, N. L., et al. (2022). Spatio-temporal level variations of the Martian seasonal north polar cap from co-registration of MOLA profiles. *Journal of Geophysical Research: Planets*, 127(10), e2021JE007158. https://doi.org/10.1029/2021JE007158

Xiao, H., Stark, A., Schmidt, F., Hao, J., Su, S., Steinbrügge, G., & Oberst, J. (2022). Spatio-temporal level variations of the Martian seasonal south polar cap from co-registration of MOLA profiles. *Journal of Geophysical Research: Planets*, 127(7), e2022JE007196. https://doi.org/10.1029/2022JE007196

Xiao, H., Stark, A., Steinbrügge, G., Hussmann, H., & Oberst, J. (2021). Processing of laser altimeter time-of-flight measurements to geodetic coordinates. *Journal of Geodesy*, 95(2), 22. https://doi.org/10.1007/s00190-020-01467-4

Xiao, H., Stark, A., Steinbrügge, G., Thor, R., Schmidt, F., & Oberst, J. (2022). Prospects for mapping temporal height variations of the seasonal CO₂ snow/ice caps at the Martian poles by co-registration of MOLA Profiles. *Planetary and Space Science*, 214C, 105446. https://doi.org/10.1016/j.pss.2022.105446

Xiao, H., Xiao, Y., Su, S., Schmidt, F., Lara, L. M., & Gutierrez, P. J. (2024). Thickness of the seasonal deposits by examining the shadow variations of the fallen ice blocks at Martian North Pole [Dataset]. *Mendeley Data*. https://doi.org/10.17632/5yy475dbry.1

Zurek, R., Campbell, B., Byrne, S., Calvin, W., Carter, L., Clancy, R., et al. (2021). Mars next orbiter science analysis Group (NEX-SAG): White paper report to the 2023-2032 planetary sciences and astrobiology decadal survey. *Bulletin of the American Astronomical Society*, 53(4), 038. https://doi.org/10.3847/25c2cfeb.da7b0e6a