Spatio‐Temporal Level Variations of the Martian Seasonal North Polar Cap From Co‐Registration of MOLA Profiles

The seasonal deposition and sublimation of CO2 constitute a major element in Martian volatile cycles. We reprocess the Mars Orbiter Laser Altimeter (MOLA) data and apply co‐registration procedures to obtain spatio‐temporal variations in levels of the Seasonal North Polar Cap (SNPC). The maximum level over the Residual North Polar Cap (RNPC) is 1.3 m, approximately half of that at the south pole (2.5 m). However, the maximum level in the dune fields at Olympia Undae can be up to 3.8 m. Furthermore, off‐season decreases up to 3 m during the northern winter at Olympia Undae are observed. These are likely due to metamorphism effects accentuated by the reduced snowfall at this period. Meanwhile, off‐season increases of up to 2 m during the northern spring are noted, the cause of which remains to be explored. The volume of the SNPC peaks at the end of northern winter and is estimated to be approximately 9.6 × 1012 m3, which is 2% more than that of the Seasonal South Polar Cap. The bulk density of the SNPC can go through phased decreases in accordance with phased accumulation at northern high‐latitudes. These findings can put important constraints on the Martian volatile cycling models.

with the cold-jetting and late re-condensation (Andrieu et al., 2018). Whereas the entire SNPC gets progressively covered by a veneer of water ice of thickness ∼200 μm during northern spring, which can partially hide the spectral signature of CO 2 ice (Appéré et al., 2011;Pommerol et al., 2013;Titus et al., 2020). It was also observed that a small amount of water ice deposition during late summer can repeatedly occur over the residual ice caps at both poles (Brown et al., 2016). Unlike at the south pole, there exist vast circumpolar dark dune fields that surround the permanent north polar cap, which are also referred to as the polar erg. Olympia Undae, which roughly spreads from 120°E to 240°E in longitude and 78°N to 83°N in latitude (refer to Figure 1 for its location), is the largest continuous dune field on Mars. Olympia Undae covers an area of approximately 4.7 × 10 5 km 2 and features linear dunes with percentage coverage typically greater than 80% (Hayward et al., 2010). The dunes at Olympia Undae can rise about 60 m above their surroundings, and feature slopes up to 35° (Putzig et al., 2014) and a horizontal decorrelation length less than 250 m (Aharonson et al., 2001). In addition, the composition of the dunes at Olympia Undae are basaltic sand overlying shallow ground ice or ice-cemented sand, leading to a low thermal inertia (Putzig et al., 2014). Through hyperspectral imaging, Olympia Undae has been identified as the largest deposit of gypsum on the surface of Mars (Das et al., 2021;Langevin et al., 2005;Massé et al., 2012). Dynamic activities associated with sublimation of the SNPC, for example, dark fans, polygonal cracks, furrows, and alcoves mostly occur within the polar erg. Although araneiform terrain is pervasive at the south pole, it has not been observed in the north Hao et al., 2020;Portyankina et al., 2013). On the steep scarps at the margins of the north polar layered deposits, mass-wasting processes like dust-ice avalanches and block falls have been observed which could be triggered by seasonal sublimation of the SNPC (Fanara et al., 2020;Russell et al., 2008;Su et al., 2021).
The direct level variation measurements of the seasonal polar caps have been made by exploiting dynamic laser altimetric profiles acquired by the Mars Orbiter Laser Altimeter (MOLA) onboard the Mars Global Surveyor (MGS) (Aharonson et al., 2004;Smith et al., 2001;Xiao, Stark, Schmidt, et al., 2022), measuring rock shadow length changes in high-resolution optical images (Mount & Titus, 2015), and Bayesian inversion involving a radiative transfer model (Andrieu et al., 2018). While the latter two approaches can be spatially and temporally limited, MOLA estimates can be extended to cover the entire polar regions. Using profile analysis at various latitudinal annuli, Smith et al. (2001) measured the maximum seasonal level variations to be ∼1 m at both poles. Aharonson et al. (2004) carried out approximations by fitting sinusoidal curves with annual and semi-annual terms to the height differences obtained at locations where two MOLA profiles intersect, that is, cross-overs, and estimated the maximum level variations to be ∼1.5 m at the north pole and ∼2.5 m at the south pole. Interestingly, Aharonson et al. (2004) demonstrated that the magnitudes of the semi-annual sinusoidal terms at Olympia Undae are unexpectedly higher (∼0.8 m) than elsewhere (∼0.2 m). Unfortunately, no additional information is available in their study to shed light on this phenomenon.
When using MOLA to measure the seasonal level variations, a common approach is to analyze the differential height measurements at cross-overs (Aharonson et al., 2004;Smith et al., 2001). However, this cross-over method can suffer from significant interpolation errors due to the large spacing between MOLA consecutive footprints (∼300 m), especially when the surface is rough with a short decorrelation length, for example, at the Olympia Undae and crater rims (Aharonson et al., 2001;Neumann, Abshire, et al., 2003). Furthermore, residual orbit, attitude, and timing errors can translate into lateral shifts of the MOLA profiles and compromise the estimation of the cross-over height differences (see also Stark et al., 2015;Steinbrügge et al., 2015;Xiao, Stark, Steinbrügge, et al., 2021). Besides, the acute angles between intersecting ground tracks, especially for laser altimeter on a polar orbit around slowly rotating bodies, pose an unfavorable geometry for the estimation of cross-over height differences (Mazarico et al., 2010;Neumann et al., 2001). For these near-parallel profiles, even small cross-track shift components are capable of shifting the intersection points dramatically upwards or downwards in the along-track The missing coverage poleward of 87°N is due to a lack of nadir-pointed profiles as a result of the orbital inclination of the spacecraft. The projection is North Polar Stereographic, with a spatial resolution of 1 km/pixel. The Residual North Polar Cap (green polygons), the polar erg dune fields including Olympia Undae (black polygons ;Skinner Jr. et al., 2006), and Chasma Boreale are marked. direction (Xiao, Stark, Steinbrügge, et al., 2021). To combat these issues, Xiao, Stark, Steinbrügge, et al. (2022) proposed the co-registration of dynamic profile segments to Digital Terrain Models (DTMs) representing a static mean reference surface as a novel approach for the detection of seasonal level variations. Meanwhile, Xiao, Stark, Steinbrügge, et al. (2022) devised a post-correction procedure based on pseudo cross-overs to correct for both a global systematic temporal bias inherited in MOLA footprint heights and all categories of residual errors. A pseudo cross-over is formed by a pair of laser profile segments that does not necessarily intersect, but instead are connected by an overlapping reference DTM (Barker et al., 2016;Xiao, Stark, Steinbrügge, et al., 2022). The height difference at a pseudo cross-over is assigned as the difference in the height corrections when the two profile segments are co-registered to the reference DTM. Based on these proposed approaches, Xiao, Stark, Schmidt, et al. (2022) temporally measured the seasonal level variations of the SSPC at grid elements of 0.5° in latitude and 10° in longitude from 60°S to 87°S. They found that maximum level can be up to 2.5 m over the Residual South Polar Cap, and its peak level increased by ∼0.5 m from Martian Year 24 (MY24) to MY25.
The masses of the seasonal polar caps have been measured by temporal gravity variations (Genova et al., 2016;Karatekin et al., 2006;Smith et al., 2009) and neutron and gamma ray flux (Feldman et al., 2003;Kelly et al., 2006;Litvak et al., 2007;Prettyman et al., 2009). Their masses have also been modeled by General Circulation Models (GCMs; e.g., Smith et al., 1999) and energy balance models (e.g., Kieffer & Titus, 2001;Schmidt et al., 2010). While consistent maximum mass estimates have been obtained for the SSPC (∼6 × 10 15 kg), the existing maximum mass estimates of the SNPC can widely range from 3.5 to 6 × 10 15 kg. By combining these mass measurements and the level variations from MOLA, the bulk density of the seasonal snow/ice can be inferred. Time-averaged bulk density for the SNPC has been constrained to be 500 ± 100 kg/m 3 for regions poleward of 85°N from Aharonson et al. (2004), 400 kg/m 3 from , and ∼600 kg/m 3 from Haberle et al. (2004). However, snow is believed to undergo densification due to self-compaction and re-crystallization, and the bulk density of the SNPC can increase from an assumed minimum of 100 kg/m 3 for fresh deposits to 1,150 ± 70 kg/m 3 at the end of northern spring (Matsuo & Heki, 2009). In contrast to the low (average) density deposits as suggested by gravity variations, neutron and gamma ray measurements, and energy balance modeling, thermal infrared and near infrared spectrometers have indicated otherwise that the SNPC is primarily composed of slab ice (Kieffer et al., 2000;Kieffer & Titus, 2001). This discrepancy is supposed to be due to a durifrost layer forming on the surface caused by basal metamorphism. This thin layer can appear as slab ice to these spectrometers which can only probe the upper centimeter of the seasonal ice caps (Kieffer, 2007;Mount & Titus, 2015).
This amount of mass being added and removed from the surface is enough to induce significant elastic displacements of the lithosphere, much like seasonal mass loading on Earth from ocean tides. Previous studies have used GCMs to characterize this load and resulting displacements (Métivier et al., 2008). Wagner et al. (2022) used results from this study and estimated up to centimeter scale deflection at the height of the northern hemisphere's winter. The spatial heterogeneity of the deflection mirrors the spatial variation of the SNPC, as expected. Characterizing this displacement is important when considering the current thermal and rheologic state of the lithosphere. Being able to detect this modeled displacement through use of current altimetry or a future geodetic mission would help further constrain the state of the Martian interior.
In this study, we focus on the north polar regions from 60°N to 87°N. By applying the proposed co-registration approach complemented with the post-correction procedure (Xiao, Stark, Schmidt, et al., 2022;Xiao, Stark, Steinbrügge, et al., 2022), we derive time-resolved seasonal level variations at grid elements of 0.5° in latitude and 10° in longitude. Characteristics of the SNPC are then extracted at individual grid elements and gridded to thematic maps for inspection. Particularly, we draw attention to some abnormal behaviors of the level variations of the SNPC within the dune fields at Olympia Undae.

Reprocessed MOLA Data Set
MOLA sampled the Martian topography, together with roughness and reflectance at 1,064 nm, at footprints of approximately 160 m in diameter with along-track spacing of around 300 m in the nadir configuration (Zuber et al., 1992). The MOLA Precision Experimental Data Record (PEDR) data set  includes more than 8,700 profiles acquired from March 1999 to June 2001, which spans from MY24 to MY25. To post-correct for residual orbit, timing, and pointing errors, Neumann et al. (2001) implemented a global cross-over analysis of the MOLA profiles with uncertainties parameterized by slowly varying harmonic functions in along-track, 10.1029/2021JE007158 4 of 21 cross-track, and radial directions. The accuracy of individual footprints has been improved to be approximately 1 m radially and within 100 m laterally ). Unfortunately, the real level variations of the seasonal polar caps were also absorbed in the aforementioned post-correction. Thus, only the original PEDR data set can be used for the purpose of this study, which was geolocated using outdated MGS trajectory and IAU2000 Mars rotational model. Therefore, we have incorporated an updated orbit model from Konopliv et al. (2006) and IAU2015 Mars rotational model as in Archinal et al. (2018) to update the MOLA geolocation (Xiao, Stark, Chen, & Oberst, 2022;Xiao, Stark, Steinbrügge, et al., 2022). Meanwhile, a pointing aberration correction, which the PEDR processing ignored, has been added during the reprocessing to account for the special relativity effects (Xiao, Stark, Steinbrügge, et al., 2021). The equation from the laser time-of-flight to footprint inertial coordinates with Mars as the observer writes, where r 1 is the one-way vector of the laser pulse outgoing leg. c is the speed of light in vacuum, τ is the laser time-of-flight. β = r 12 /τc, in which r 12 denotes the vector directed from the spacecraft's position at the laser transmission time t 1 to that at reception time t 1 + τ. e 1 is the normalized boresight vector denoting the orientation of the emitted laser pulse. r in is the positional vector of the bounce point with respect to the Center-of-Mass (CoM) of Mars in the International Celestial Reference Frame (ICRF). r s (t 1 ) is the spacecraft's positional vector with respect to the CoM of Mars at t 1 . The obtained inertial coordinates of the footprint r in can be further converted to body-fixed coordinates by evaluating the Mars rotational model at laser reflection time t 1 + r 1 /c (Xiao, Stark, Steinbrügge, et al., 2021). For the nadir-pointing profiles, the pointing aberration can result in lateral shifts of up to 5 m and radial ones up to 5 cm in magnitude (Xiao, Stark, Steinbrügge, et al., 2021). Due to the opposite sign in the radial direction for ascending and descending tracks, this factor can actually lead to height variation offsets of magnitudes up to 10 cm. The self-consistency of the reprocessed MOLA PEDR data set, evaluated as the root-mean-square height discrepancy at cross-overs in equatorial and tropical regions where features no seasonal CO 2 deposition, is significantly better than the original one (Xiao, Stark, Steinbrügge, et al., 2022). The root-mean-square value has been reduced from 8.2 to 7.0 m, a nearly 15% improvement. But, there still remain significant profile shifts, including a global systematic temporal bias in height (see also Section 3; Xiao, Stark, Steinbrügge, et al., 2022), due to timing, pointing and other unidentified errors. These error sources need to be taken care of in order to accurately and precisely obtain the seasonal level variations at the polar regions (Xiao, Stark, Chen, & Oberst, 2022;Xiao, Stark, Steinbrügge, et al., 2022).

MOLA Reference DTM From Self-Registration
To build a reliable reference DTM of the Martian mean surface at [40°N, 87°N] as preparation for the extraction of the seasonal level variations, we resort to self-registration of the reprocessed MOLA profiles and a divide-andconquer strategy as in Xiao, Stark, Schmidt, et al. (2022). When the lengths of the profiles are too short, the co-registration can be unstable. Meanwhile, when they are too long, there exist significant approximation errors, as we are assuming constant corrections in line, sample, and height, that is, map coordinates, in polar stereographic projection. For this reason, the self-registration is separately applied in five annular regions, namely, [ ]. The self-registration process is iterated until the refinement is visually unidentifiable. First, we randomly select 25% of the reprocessed MOLA profiles. Then, we generate an intermediate DTM with a spatial resolution of 1 km/pixel from the rest of the profiles (75% of the total). After that, we individually co-register the selected profiles to the intermediate DTM with a parameterization in map coordinates. We iterate the aforementioned process 30 times until all outlier profiles are removed and speckle noise cannot be further reduced. Finally, a MOLA reference DTM with a spatial resolution of 1 km/pixel can be obtained by gridding these self-registered altimetric profiles.

Mapping of Seasonal Level Variations of the SNPC
In this study, we use a method proposed by Xiao, Stark, Schmidt, et al. (2022) and Xiao, Stark, Steinbrügge, et al. (2022) to map the seasonal level variations by taking advantage of height corrections from local co-registration processes. The seasonal level variations of the SNPC are extracted in a spatial grid of elements of 0.5° in latitude and 10° in longitude. First, a grid element, a laser profile within this grid element, and one laser footprint along the track are chosen. Then, all footprints acquired within ±30 s of the acquisition time of the footprint in question are selected to form a local profile segment. For nadir-pointing profiles, the selected profile segment contains about 300 footprints on each side and can extend to about 180 km. For comparison, the profile within each grid element is on average ∼90 footprints in length. Then, the selected profile segment is co-registered to the MOLA reference DTM parameterized in stereographic map coordinates (Xiao, Stark, Steinbrügge, et al., 2022). The height correction from the aforementioned local co-registration is considered as the observable, which stands as the height variation with respect to the reference DTM. Time stamp associated with this observable is assigned as the acquisition time of the central footprint. The aforementioned process is then repeated for all footprints that fall within the grid element. Finally, all observables are median-binned against their time stamps with a temporal resolution of ∼5 days, to form a level variation time series. The precision is represented by the scaled median absolute deviation from the median (Leys et al., 2013). Through the co-registration process, the disadvantages of using the height differences at cross-overs, that is, errors due to interpolation and the misalignment between the profiles, can be alleviated. Exceptionally, for grid elements situated within Olympia Undae, we extend the length of the local profile segments from 180 to 360 km to stabilize the co-registration processes. The roughness of these dune fields at Olympia Undae is high at small scale, for example, within the MOLA footprint (Neumann, Abshire, et al., 2003), but extremely low at large scale (Kreslavsky & Head III, 2000). Without large-scale features, the co-registration can fail and yield erratic results.

Post-Correction Procedure
Meanwhile, there exists a global systematic temporal height bias in the MOLA geolocation even after its reprocessing. This temporal bias does not vary much with latitude and features a peak-to-peak amplitude of 2 m with a phase that resembles the synodic period of Mars (Xiao, Stark, Steinbrügge, et al., 2022). To correct for this artifact while aiming for high-precision temporal level variations, we exploit a process known as regional pseudo cross-over adjustment (RPCA) which combines the concept of DTM-based pseudo cross-over and the local co-registration technique (Xiao, Stark, Steinbrügge, et al., 2022). For the pseudo cross-over formed by the kth and lth profile segments, the height misfit at the pseudo cross-over (Δh pc ) is assigned as where Δh k and Δh l are the height corrections from the co-registration of the kth and lth profile segments to the reference DTM with a parameterization in polar stereographic map coordinates, respectively. The RPCA is similar to the regional cross-over adjustment (RCA) that is carried out over a relatively small, non-dynamic area to improve the internal consistency of the biased altimetric profiles on Earth (Ewert et al., 2012). However, the RPCA is advantageous over the RCA in multiple aspects. As any combination of a pair of profiles can form a pseudo cross-over, the available number of "cross-overs" is normally several times that of the traditional cross-overs. The co-registration processes in the formation of the pseudo cross-overs can avoid the interpolation errors and compensate for lateral shifts of the profiles associated with the usage of traditional cross-overs. The unfavorable shallow intersecting angles when calculating the height misfits at the traditional cross-overs can be avoided. Furthermore, the profile segment pairs that constitute pseudo cross-overs can be far apart across the research region, for example, a latitudinal annulus, to impose as "global" constraints. In contrast, the profiles that intersect to form the traditional cross-overs have to be spatially proximal. We parameterize a constant height correction for each of the profile segments and solve for the unknowns by exploiting a least squares that minimizes the root-mean-square height misfit at pseudo cross-overs. Ridge regression is adopted in the RPCA to solve the ill-posed least squares problem due to the multicollinearity between the large number of model parameters (Hoerl & Kennard, 1970;Xiao, Stark, Steinbrügge, et al., 2022). In general, the regularization provides improved efficiency in parameter estimations in exchange for a tolerable amount of biases, that is, bias-variance trade-off.
In practice, we consecutively apply two RPCAs in a post-correction procedure as "bi-RPCA" to associate the polar regions with CO 2 level variation signal with a tropical annulus of [44°N, 56°N] that is assumed to be free from seasonal level variations (Xiao, Stark, Steinbrügge, et al., 2022). First, each of the profile segments intercepted by annulus [44°N, 56°N] is co-registered in stereographic map coordinates to the reference DTM to obtain height corrections per profile for the global systematic temporal bias. The first RPCA is then applied to acquire the adjustments for the aforementioned height corrections. Assuming the systematic temporal bias is spatially invariable and equals that acquired at [44°N, 56°N], these height corrections are then adjusted and applied to correct the profile segments at the polar regions. Subsequently, the second RPCA is performed at a grid element at the north pole to obtain the height corrections for each of the profiles due to all categories of residual errors, including that of the reprocessed MOLA profiles, the co-registration process, and the reference DTM. The local profile segments used in the second RPCA process are formed by selecting footprints acquired within ±30 s of the central footprints along the profiles intercepted by the grid element. Finally, height corrections from both of these RPCAs are applied to observables from Section 3.2, which are then median-binned to form a post-corrected seasonal level variation time series. It should be noted that while pseudo cross-overs acquired within 5 days are used in the first RPCA, those acquired within 15 days are adopted in the second RPCA to boost their numbers so as to further suppress the noise. The maximum level variation rate at the north pole is generally much smaller than that at the south pole (Aharonson et al., 2004), so real level variations that occur within 15 days do not significantly interfere with the post-correction procedure.

Gridded Products
We apply the proposed techniques to all grid elements from 60°N to 87°N. Key features, that is, the peak-to-peak level and early height accumulation at around L s 300° (marked in Figure 4), are then extracted from the obtained temporal level variation in each element and gridded to thematic maps (available at Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2021)). For grid elements from 78°N to 84°N and 110°E to 250°E that encompass Olympia Undae, we additionally examine the spatial distribution of the off-season decreases in northern winter and the off-season increases in northern spring (marked in Figure 5). To suppress noise, these gridded maps are filtered using a local sliding 3 × 3 Gaussian window (1.5° × 30°) with a sigma of five grid elements in both longitudinal and latitudinal directions. Smooth interpolation using spherical splines with global gradient estimates, and without tension, is also applied to assign values within the grid elements.

Results
The reference DTM at the Martian north pole produced from the self-registered MOLA profiles is shown in Figure 1 (available at Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2021)). Locations of mentioned geological units and topographic features are marked. The reference model covers areas poleward of 60°N at a spatial resolution of 1 km. It is used in the local co-registration processes in deriving the temporal level variations and the corresponding post-corrections.
To obtain the general trend of the temporal level variations as a function of latitude, the local co-registration approach and the bi-RPCA procedure are also applied to seven different annuli (centered at 48°N, 62°N, 72°N, 80°N, 82°N, 84°N, and 86°N with a latitude width of 2°). Pseudo cross-overs acquired within 5 days are adopted in both of the RPCAs in the bi-RPCA procedure. The comparison between the results is shown in Figure 2. The precision of these results can be indicated by the degree of scattering of the measurements, which is on the level of several centimeters. The accumulation begins at early to middle northern fall, much more advanced tha at the south where it normally begins in southern winter (Figure 2a in Xiao, Stark, Schmidt, et al. (2022)). For temporal level variation curves along annuli 80°N and 82°N, the maximum levels are reached at around the end of the northern fall. For the other annuli, the levels peak at the end of northern winter, which is much later than that for annuli 80°N and 82°N, after which they begin to decrease at northern spring due to sublimation. The magnitude of the maximum level variations gradually decreases from ∼1.4 m at 86°N to ∼0.8 m at 72°N, and then to ∼0.2 m 7 of 21 at 48°N. As the lowest latitude the SNPC can extend to is approximately 50°N, ∼0.2 m at 48°N can be considered to some extent as the accuracy of the obtained temporal level variations at these annuli. Two previous results from MOLA measurements, that is, Smith et al. (2001) at latitudinal annulus 86.5°N (green line) and Aharonson et al. (2004) at latitudinal annulus 86°N (light blue line) are shown for comparison ( Figure 2). Our curve at annulus 86°N is generally consistent with these two previous studies. However, our result at annulus 86°N does not feature significant off-season increases in the northern summers as observed by Aharonson et al. (2004). The off-season decrease in the northern fall as stated to be the result of atmospheric warming due to regional dust storms by Smith et al. (2001), is not revealed in our results. Interestingly, the level at annulus 86°N peaks at early northern winter in Aharonson et al. (2004) while that comes much later as indicated by Smith et al. (2001) and our study. Possible explanations of these differences are discussed in Section 5.1.
Further, the annulus analysis is extended to continuously cover the latitudes from 50°N to 87°N with an interval of 2° in latitude. For the gap poleward of 87°N, we assume that the snow/ice variations are identical to that acquired at the neighboring annulus 86°N. Volume for each annulus is then calculated by multiplying the snow/ice level by the area of that annulus. The latter is evaluated by retaining eccentricity terms up to degree 10 with Mars assumed as an ellipsoid with an equatorial radius of 3,396.19 km and a north polar radius of 3,373.19 km (Text S1 in Supporting Information S1). The total volume evolution of the NSPC at each time step is shown as the blue dots in Figure 3. To take care of the data gaps due to solar conjunctions, we fit the volume time series using the Gaussian Process regression in which a radial basis function kernel is adopted. The sigma added to the diagonal of the kernel matrix is assigned as three folds of its associated scaled median absolute deviation (refer to the error bars for the blue dots). The fitted curve is represented by the blue curve and 95% confidence interval is represented by the blue shade. The volume of the SNPC begins accumulating at the beginning of fall and peaks at the end of northern winter with 9.6 × 10 12 m 3 , which is at approximately the same level as that of SSPC in MY24 (9.4 × 10 12 m 3 ; Xiao, Stark, Schmidt, et al., 2022). Then, the snow/ice gradually sublimates away throughout the northern spring and summer. The volume at the middle of the northern summer in MY25 is close to 4 × 10 12 m 3 which doubles that in MY24. This abnormality is due to the artifact observations of significant snow/ice level persisting to the northern summer in MY25, as discussed in Section 5.2.
Movie S1 in Supporting Information S1 on the spatio-temporal level variations of the SNPC from 60°N to 87°N is available as the supplementary file (also downloadable at Xiao, Stark, Schmidt, Hao, Steinbrügge, et al. (2021)). We can see that after the areal extent peaks at the end of fall and begins to decline thereafter (Piqueux et al., 2015), the accumulation is still happening at high latitudes throughout the polar night in the northern winter (refer to a brief discussion in Section 5.5). A typical level variation time series for a grid element outside of the dune fields at Olympia Undae ([85.75°N, 86.25°N, 30°E, 40°E]) is illustrated in Figure 4. The level starts to increase at the beginning of northern fall and an early surface accumulation of ∼0.8 m thick at around L s 300° in the beginning of northern winter can be observed. The time stamp of L s 300° is chosen for comparison purposes with that at the south pole. The level maximizes at ∼1.3 m at the end of the northern winter. The off-season increases in the northern summers observed in Aharonson et al. (2004), that is, the light blue curve in Figure 4, are not present in our result for this specific grid element and that of Smith et al. (2001). It is worthy to note that an ice layer of ∼0.5 m in level persists all the way into late northern summer in MY25, while no significant ice layer is observed in the contemporary period in MY24. This unexpected difference in the summer level is discussed in Section 5.2.
The temporal level variations in grid elements in the dune fields at Olympia Undae are distinct from that at the other regions. To clarify, a typical level variation time series for a grid element within Olympia Undae ([81.25°N, 81.75°N, 180°E, 190°E]) is shown in Figure 5. Up to 0.5-1 m deep of snow/ice is observed in the northern summers in both MY24 and MY25. These are considered as artifacts as discussed in Section 5.2. In comparison, a limited level of snow/ice remains in the corresponding periods for the examined grid elements outside of Olympia Undae. The peak-to-peak variation is approximately 2.8 m, which is much larger than ∼1.3 m outside of Olympia Undae (Figure 4). This is unexpected and is discussed in Section 5.3. The most prominent feature is its M-shaped central portion. Instead of continuously growing, the level unexpectedly drops in the early to-middle northern winter after peaking at the end of northern fall, and then recovers during late northern winter and the beginning of northern spring. These off-season decreases at northern winter and increases mainly at northern spring, as marked in Figure 4, are elaborately examined and discussed in Section 5.4.
The map of the peak-to-peak level variation, examined in the southern winters of MY24 and MY25, is shown in Figure 6a with the associated precision map in Figure S1 in Supporting Information S1. These maps clearly illustrate the locality of the maximum level variation and its non-axisymmetric pattern around the north pole. In particular, maximum level is higher in lowlands southward of Chasma Boreale than that at the same latitudes. The largest magnitudes up to ∼3.8 m are located in the linear dune fields at Olympia Undae (top panel in Figure 6a). Elsewhere, magnitudes up to ∼1.3 m can be spotted in some patches over the Residual North Polar Cap. The precision as measured by the scaled median absolute deviation is normally better than 15 cm except at Olympia Undae where it can degrade to up to 40 cm. The peak-to-peak level variation map for SSPC from Xiao, Stark, Schmidt, et al. (2022) is also shown for comparison (Figure 6b), it can be seen that maximum levels at the north pole are generally lower than that at the south pole, except at Olympia Undae. In addition, the patch with maximum level more than 1 m in magnitude is much more elongated at the north pole than at the south pole. Statistics of the maximum snow/ice level, that is, mean, maximum, and standard deviation, are examined in each annulus against latitude as shown in Figure 7. The counterparts for SSPC from Xiao, Stark, Schmidt, et al. (2022) are also plotted as dots for comparison. The most prominent features are the bumps from 77.5°N to 83.5°N for all of the three statistics. These bumps are due to the significantly larger level variations in the dune fields within Olympia Undae as compared with elsewhere at the north pole. The abnormal behavior of the SNPC observed at Olympia Undae is thoroughly discussed in Sections 5.3 and 5.4. Evolution of the maximum peak-topeak level, which a crest centered around the latitudes of the Olympia Undae, is consistent with the map product shown in Figure 6a. Unlike that at the south pole, the mean level does not exhibit quasi-linear patterns. The curve runs almost flat at ∼0.5 m from 60°N to 67°N, which could suggest that the level variations less than 0.5 m in magnitude are concealed by the residual biases in the MOLA measurements. Thus, 0.5 m can be the approximate accuracy of the temporal level variations obtained at individual grid elements. The variability of the snow/ice level within each latitudinal annulus can be indicated by the standard deviation, which stays within 15 cm except at latitudes that encompass portions of the dune fields at Olympia Undae.
The early height accumulation map measured at around L s 300° in northern winter of MY24 at the north pole is shown in Figure 8a and Figure S2 in Supporting Information S1. It can be observed that the radial expansion of the SNPC is spatially heterogeneous. Magnitudes greater than 1 m, that is, the reddish patches, cluster at the regions surrounding Chasma Boreale and within Olympia Undae. However, it is worthy to note that the maximum level at Olympia Undae is actually reached at around the northern winter solstice, and off-season decreases have already taken place before L s 300° ( Figure 5). The blueish-to-greenish transition outline, corresponding to early height accumulation of ∼0.7 m, roughly coincides with the exterior borders of the polar erg ( Figure 1). The early height accumulation map measured at around L s 120° in southern winter of MY24 at the south pole, adapted from Xiao, Stark, Schmidt, et al. (2022), is also shown for comparison (Figure 8b and Figure S2 in Supporting Information S1). It can be seen that relatively . Each black dot in the background represents a bi-regional pseudo cross-over adjustment post-corrected height correction from the local co-registration of a profile segment centered at a specific footprint. Yellow shades denote the precision of the temporally binned median temporal curve measured by the scaled median absolute deviation. Two key statistical characteristics for interpretation purposes, that is, maximum level of the seasonal polar cap and early winter level at around L s 300°, are marked, which are extracted at all individual grid elements and mapped to be spatially examined. Two previous results of Smith et al. (2001) at annulus 86.5°N and Aharonson et al. (2004) at annulus 86°N with precision represented by shades are also shown for reference. large magnitudes cluster in limited regions. Overall, early height accumulation at the north pole is higher and more concentrated as a whole than that at the south pole. These observations contrast the fact that the maximum level of the SNPC is in general significantly lower than that of the SSPC ( Figure 6). This contradiction could be partially attributed to a longer and colder southern winter due to the eccentricity of the Martian orbit (e = 0.09) and the near-coincidence of its apocenter and southern winter solstice.

Comparison to Previous Results at Annuli 86°S/N
The availability of temporal level variations at annuli 86°S/N from Smith et al. (2001) and Aharonson et al. (2004) enables a direction comparison to our results (Figure 2 for annulus 86°N and Figure 2 in Xiao, Stark, Schmidt, et al. (2022) for annulus 86°S). While the maximum level at annulus 86°N from Aharonson et al. (2004) is similar to ours, their magnitude at annulus 86°S is significantly larger than ours. Previously, we have demonstrated that both reprocessing of the MOLA profiles and the RPCA processes do not lead to systematic temporal biases to the footprint heights (Xiao, Stark, Steinbrügge, et al., 2022). Thus, the difference between our results and that of Aharonson et al. (2004) at annuli 86°S/N could mainly be attributed to different post-correction strategies utilized to correct for the global systematic temporal bias (Xiao, Stark, Steinbrügge, et al., 2022). We have done the correction assuming the systematic temporal bias to be constant at the polar regions and equals that obtained at annuli 50°S/N. In contrast, Aharonson  Thus, it is imperative that we pin down the root cause for this global temporal bias and correct the MOLA heights for it on a local basis. Smith et al. (2001) differentiated the level variation signal in each polar annulus to that at 60°S/N to correct for the global temporal bias, similar to our approach. Although the peak-to-peak level variation at annuli 86°N from Smith et al. (2001) is comparable to ours and that of Aharonson et al. (2004), their result is surprisingly lower in magnitude at annulus 86°S. This may arise from additional temporal smoothing that they have applied as their curves are the smoothest among the three results compared here.

Summer Level in MY25
From the level variation time series for the grid element of [85.75°N, 86.25°N, 30°E, 40°E] (Figure 4), we can see that a ∼0.5 m snow/ice level remains to the northern summer of MY25. The distributions of seasonal snow/ice level at L s 110° in the northern summers of MY24 and MY25 are shown in Figure 9 and Figure S3 in Supporting Information S1. In contrast to that in MY24, significant portions within 70°N to 80°N feature ∼0.5 m thick of snow/ice deposits in the northern summer of MY25. An exception is in the dune fields at Olympia Undae where the SNPC features rather high summer level (0.5-1.5 m) in both MY24 and MY25. Wang and Ingersoll (2002)   made global maps of the north polar regions of Mars using the MGS Mars Orbiter Camera (MOC) wide-angle swaths during the northern spring and beginning of northern summer in MY25. Unfortunately, no apparent seasonal deposits can be observed in these maps outside of the Residual North Polar Cap at L s 111°. In addition, no seasonal snow/ice deposits at latitudes between 70°N and 80°N were observed by the Mars Color Imager (MARCI) onboard Mars Reconnaissance Orbiter (MRO) during the northern summers at MY29, MY30, and MY31 (Calvin et al., 2015). Thus, these MOLA-measured 0.5-1.5 m level of summer snow/ice is considered as an artifact, which contributes to the non-diminishing trend of the SNPC volume from L s 30° in the northern spring to the end of the northern summer in MY25 (Figure 3).
In the dune fields at Olympia Undae, there indeed exist bright patches in between the dunes that can linger throughout the northern summer ( Figure 12 in Section 5.4). Majority of these white patches cannot be effectively resolved due to the coarse resolution of these global MOC maps (7.5 km/pixel) and MARCI images (∼1.9 km/ pixel). As discussed later in Section 5.4, these bright patches are most likely made of gypsum and limited evaporitic minerals that can persist throughout the entire northern summer. Thus, the up to 1.5 m level of snow/ice observed at Olympia Undae during the beginning of northern summers in both MY24 and MY25 is also considered as an artifact.
Unfortunately, we do not have direct evidence to prove the signals extracted at seasons other than the northern summers are correct. A supporting evidence is that the level obtained during northern fall, winter, and spring orderly and gradually decreases from ∼1.3 m at annulus 86°N to nearly 0 at annulus 48°N ( Figure 2). Exceptions are the curves for annuli centered at 80°N and 82°N, where the accumulation comes more rapidly at mid-tolate northern fall. However, this can be expected as these annuli encompass portions of Olympia Undae, where the level of the SNPC peaks early at the end of the northern fall ( Figure 5). The close to zero level, although with limited deviations of ∼0.2 m due to residual errors, at annulus 48°N is consistent with the observations that the maximum SNPC can roughly extend down to 50°N (Piqueux et al., 2015). In addition, we would like to mention that low-elevation areas adjacent to Chasma Boreale feature higher level than regions at the same latitudes between 70°N and 80°N (Figure 6a). This phenomenon is consistent with the neutron spectroscopy observations that CO 2 deposition is skewed toward areas exposed to frigid katabatic winds from Chasma Boreale (Prettyman et al., 2009). Interestingly, we have also observed systematic non-zero levels during southern summer and fall for the SSPC at the south pole (Xiao, Stark, Schmidt, et al., 2022). For that, we have presented several strong arguments why these off-season biases do not act on the seasonal signals extracted during the southern winter and spring. This serves as another evidence that the artifacts observed in the northern summers should not significantly undermine the validity of level variations at the other seasons.

High Maximum Snow/Ice Level at Olympia Undae
The peak-to-peak level at Olympia Undae dune fields is much larger than at the other regions at the north pole ( Figure 6). The composition of the dunes at Olympia Undae is basaltic sand overlying shallow ground ice or ice-cemented sand, resulting in low thermal inertia (Putzig et al., 2014). Thus, the condensation rate at northern fall and beginning winter, and hence the maximum snow/ice level, can be much higher at these dune fields. The low thermal inertia of the dune fields combined with their vast expansion with high dune percentage can create a local climate that favors the surface condensation and atmospheric deposition than the other dune fields within the polar erg. The high rise (up to ∼60 m) of the dunes at Olympia Undae and their high latitudes can also enable extensive shadowing at pole-facing slopes and in between the dunes during the deposition seasons, facilitating the CO 2 frost condensation (e.g., Khuller et al., 2021). Calvin et al. (2015) made an animation of the retreat of the SNPC during the northern spring and summer of MY29, MY30, and MY31 using temporal MARCI images. Interestingly, the complete removal of the white patch at Olympia Undae significantly lagged behind that in regions at similar latitudes (see also Figure 18 in Cantor et al., 2010). This phenomenon stands as a supporting evidence of a much higher snow/ice level at Olympia Undae at the end of northern spring. One thing noteworthy in the MARCI animation is that the last patch at Olympia Undae to survive before being sublimated away lay at ∼210° in longitude, which is spatially consistent with the obtained maximum level map ( Figure 6). However, no exceptionally smaller emissivity, that is, greater quantities of snowfall, within Olympia Undae was observed compared to other regions during the northern fall from MY29 to MY33 (Gary-Bicas et al., 2020). This contradiction could possibly be due to higher direct condensation at Olympia Undae or inter-annual difference in the quantities of snowfall as the MOLA measurements were acquired in MY24 and MY25.

Off-Season Decreases and Increases at Olympia Undae
Dune fields within Olympia Undae typically exhibit M-shaped temporal level variation curves ( Figure 5). The off-season decreases during early to-middle northern winter and the subsequent off-season increases during late northern winter and beginning of northern spring can be up to three m ( Figure 10 and Figure S4 in Supporting Information S1). One aspect to consider is the saturation-related issue of MOLA due to the absence of a variable gain amplifier and limited digital ranges in width and amplitude recordings Xiao, Stark, Schmidt, et al., 2022). Saturation can bias the level variation curve upwards with a median of 0.64 m due to the minimum estimates of the leading edge correction scheme . Fortunately, due to the dark materials of the dune fields, majority of the return pulses were not saturated at Olympia Undae. This unsaturation can be illustrated in a comparison of the return pulse energies within a grid element over the Residual North Polar Cap ([85.75°N, 86.25°N, 30°E, 40°E]) and a grid element at Olympia Undae ([81.25°N, 81.75°N, 180°E, 190°E]) ( Figure 11). Emitted pulse energy mainly varied with environmental temperature, declining overall with age. Meanwhile, frequent rapid fluctuations can also be observed which could be due to failures of individual laser diode bars . For the grid element over the Residual North Polar Cap, the received pulses are saturated most of the time except during periods in northern fall and northern winter, when phased CO 2 snowfalls are expected which can scatter and attenuate the laser pulses (Gary-Bicas et al., 2020). These phased snowfalls inferred from cloud opacity and emissivity from MCS measurements in MY29 to MY33 correspond well with the absorptive clouds as revealed by MOLA measurements in MY24 and MY25 . In contrast, a majority of the received pulses are unsaturated for the gird element over dark dune fields at Olympia Undae. In light of this, we also solve for the product of reflectivity (assuming Lambertian) r s and two-way Mars atmosphere transmissivity 2 t , or normalized "return energy", using MOLA's link equation (blue dots in Figure 11; Abshire et al., 2000) where E r is the echo pulse energy, E tr is the transmitted laser-pulse energy, τ r is the receiver optics transmission, A r is the receiver telescope aperture area, and R m is the laser one-way range. The product denotes return strength with respect to that from a perfectly reflecting Lambertian surface in clear atmosphere. The maximum of this normalized "return energy" (∼0.3) is reached at around L s 60° in MY25 which is followed by a sharp decrease till the end of the northern spring. The decrease of the (normalized) return energy during the early to-middle northern winter can be due to the increasing translucency of the snow/ice layer as a result of the metamorphism effects. In addition, the increasing trend during late northern winter and early northern spring can be due to the formation of an overlying water ice veneer. Pommerol et al. (2013) utilized the Compact Reconnaissance Imaging Spectrometers for Mars (CRISM) to examine the CO 2 and H 2 O ice band depth at Zanovar (81°N, 156°E) within Olymipa Undae (also compare to Figure 26 in Hansen et al. (2013) of the temporal evolution of this dune field observed using High Resolution Imaging Science Experiment images (HiRISE; McEwen, 2005;McEwen et al., 2007). The CO 2 ice band depth has been shown to continuously decrease and the H 2 O ice band depth continuously increase during the northern spring (data coverage from approximately L s 20° to L s 60° in MY29 and MY30). For unsaturated pulses, the uncertainty for the range walk corrections, using half of the measured pulse width, is typically ∼30 cm and these random errors can be largely compensated for during the large averaging in the co-registration processes (up to 1,200 footprints in each profile segment) and the RPCA post-correction procedure. Thus, saturated-related biases are unlikely to be responsible for these off-season decreases and increases observed within Olympia Undae. Another aspect to consider is the global systematic temporal bias, which can spatially fluctuate but treated as a constant that equals the one obtained at latitudinal annulus of [44°N, 56°N] (Figure 25 in Xiao, Stark, Steinbrügge, et al. (2022)). Due to its similarity to the synodic period of Mars, the systematic temporal bias is mostly likely a result of untreated biases in MGS orbit, attitude, or thermally induced MOLA alignment, and should not be affected by the specific surface properties, for example, sand dunes at Olympia Undae. Besides, the magnitudes of its fluctuations with latitudes are normally below 0.5 m, not comparable to the off-season decreases and increases that can be up to 3 m.
At around L s 160°-165° in late northern summer, thin haze rapidly develops into the thick water ice clouds known as the north polar hood. The polar hood can last from late summer, fall, and all the way to winter and early spring (Benson et al., 2011;Brown et al., 2016;Calvin et al., 2015;Navarro et al., 2014). The distribution of this optically thick polar hood encompasses all ranges of longitudes and does not concentrate over Olympia Undae (refer to Figure 6 of Wang & Ingersoll, 2002, for its distribution in MY24 and MY25). Thus, if these anomalies at Olympia Undae are due to laser pulses being reflected or attenuated by the polar hood, then similar patterns should also be revealed at other polar regions. In contrast, we find those off-season anomalies are spatially confined to Olympia Undae. Additionally, there also exists the south polar hood and we have not observed these kind of anomalies that can be up to 3 m at the south pole (Xiao, Stark, Schmidt, et al., 2022). To examine possible diurnal effects of accumulation and formation of the water ice clouds, we also check the temporal evolution of the local time at the MOLA measurements. The curves for both ascending and descending tracks turn out to be smooth-changing with variations of less than 2 hr, which can be attributed to MGS being placed in a frozen sun-synchronous quasi-circular orbit. The diurnal effects can thus be ruled out. In fact, during selection . For the latter grid element, the normalized "return energy" for the return pulses is also calculated, as shown in blue dots.
of valid ground returns, we have already used the shot quality descriptor flag to filter out returns from the clouds . In summary, the interference of the polar hood with the laser pulses should not be responsible for the off-season anomalies at Olympia Undae.
We also have to consider the high sand-flux, with dune migration rate up to ∼2 m/MY, at Olympia Undae as driven by summer katabatic winds modulated by the receding SNPC, large topographic relief, and large thermal contrast (Bridges et al., 2012;Chojnacki et al., 2019). The dune migration directions are determined by the katabatic winds and steered zonally by Coriolis-force winds, which can spatially vary . Sand dunes can become ice-cemented throughout northern fall, winter, and spring, which inhibits their migration. Nevertheless, we introduce various lateral shifts with a rate of 2 m/MY to the MOLA profiles to simulate the impact of the dune migration. We apply the proposed techniques with co-registration implemented merely in the radial direction so as not to absorb the imposed lateral shifts. However, the general trend of the obtained level variation curves stay largely unchanged and the off-season abnormalities remain. The non-repeatability of the MOLA tracks, that not re-observing the same surface periodically, and their non-alignment with the running direction of the dunes mean that the height biases induced by the dune movement can have different signs at stoss slopes and slip-faces. Then, these biases can largely average out along extending profile segments (1,200 footprints and 360 km in length) used in the co-registration processes. Thus, aeolian-driven dune migration is unlikely to account for these observed off-season decreases and increases. Another interesting form of sand mass flux is the collapses of the dune slip-faces, or dune alcoves. However, we also find this phenomenon to be unlikely cause for the anomalies observed within Olympia Undae: The alcoves are formed between early northern fall and early northern spring (Diniega et al., 2019;Hansen et al., 2015), which are probably due to frost-related activities during deposition and sublimation of the seasonal ice cover. Thus, the same effects of the alcove formations should be present throughout northern fall to spring, which are unlikely to create these vibrating M-shaped level variation curves observed at Olympia Undae; like the dune migration, these sand collapses do not change the average height, and hence total volume, of the sand dunes. As MOLA heights essentially represent averages within the laser-illuminated footprints, this mere redistribution of sand is unlikely to be responsible for these systematic variations in surface heights over the enormous Olympia Undae. Instead, dedicated investigations and realistic simulations to quantify how the MOLA pulses and resultant height measurements would be affected by varying surface due to, that is, migrating dunes, polygonal cracks, furrows, alcoves, and avalanches, in these dynamic and expansive dune fields would be needed in the future. Interaction with the varying albedo, slope, and roughness of these dynamic dunes can lead to distorted non-symmetric return pulses of, for example, broadening trailing edge, violent variations, and multiple peaks. In contrast, the MOLA range walk correction scheme assumed the reflected laser pulses to be perfectly Gaussian-shaped and the pulse width (halved) used for the purposes of correcting the unsaturated footprints were evaluated at crossings with prefixed energy thresholds. Thus, there should exist model-dependent biases in the applied range walk corrections, and hence in the MOLA heights. A viable way to check this hypothesis would be to take footprint-scale topography measured from high-resolution stereo-photoclinometric DTMs (e.g., HiRISE DTMs; McEwen et al., 2007), synthetically add topographic and albedo variations due to dynamics in the dune fields, simulate the full waveform of the returned pulses (e.g., Budge et al., 2006), and finally quantify the biases in the MOLA range walk corrections. In this way we can confidently conclude if these anomalies observed within Olympia Undae are due to artifacts in the MOLA data set.
As an attempt to verify these off-season decreases and increases, we examine a series of temporal HiRISE red band images of a typical dune field located at 81.6°N and 178.8°E during northern spring and summer of MY32 ( Figure 12). Unfortunately, no HiRISE images are available at this site in the dark northern winter. The accentuated topography at L s 1° in the beginning of northern spring indicates the presence of slab ice which forms via radiative cooling and thus conforms to the underlying topography. Then, at L s 35°, dark dust fans and polygonal cracks begin to form along the crestlines and sun-facing slip-faces where higher solar insolation is available than elsewhere. Further into the northern spring at L s 62°, these features become more pervasive and the dark erupted materials are carried away by winds to cover extensive portions of the snow/ice surface. These spring activities in the northern dunes have been previously demonstrated to be consistent with the Kieffer's model which was initially proposed for the southern hemisphere Portyankina et al., 2013).
Then at the beginning of northern summer at L s 91°, the crests of the dunes are defrosted and the erupted substrate materials are blended back to the source materials. Only a thin seasonal ice deposit layer remains as observed at L s 91°. Pommerol et al. (2013) also examined the CRISM cubes at Zanovar (81°N, 156°E) at L s 89.1° shortly before northern summer solstice, strong spectral signatures of water ice are detected as a large amount of them still persist on the north-facing flanks of the dunes (Figure 10 in Pommerol et al. (2013)). It is also observed at Olympia Undae that the CO 2 ice spectral signature is gone after L s 80° and the H 2 O ice spectral signature dominates (Appéré et al., 2011). The sources for these water ice can be sublimation of the water ice annulus at the edge of the seasonal polar cap, icy regolith, and that released during the sublimation of the seasonal ice cap (Appéré et al., 2011). Thus, the deposits observable at L s 91° could likely be a thin layer of water ice.
Subsequently at L s 114°, the aforementioned thin water ice layer also sublimates away and reveals lingering patches of bright materials between the dunes, which can persist all the way into the late northern summer (see that at L s 174°). Currently, there is no dedicated research as to the composition of these bright patches, but some early investigations suggested that they show more gypsum than water ice signatures (Fishbaugh et al., 2007;Horgan et al., 2009). But the fact that these bright patches can persist throughout the year brings another evidence that they must be some kind of mineral. In addition to gypsum, there might be other evaporitic minerals as small interannual variations of these bright patches have also been observed in the HiRISE images. These bright materials do not exhibit signs of significant areal shrinkage in the northern summer and that they cover just a fraction of the dune surface at Olympia Undae. Thus, no significant level variations from these materials are expected during the northern summer. We can then infer that those up to 0.5-1 m level observable in Figure 5 is most likely an artifact, as also discussed in Section 5.2.
One plausible hypothesis for the off-season decreases in the early to-middle northern winter is due to snow/ ice metamorphism that includes both the compaction by its own gravity and re-crystallization (Eluszkiewicz et al., 2005;Matsuo & Heki, 2009;Mount & Titus, 2015). As the build-up of the snow/ice happens rapidly during the northern fall prior to the northern winter, the self-compaction can arrive early. Additionally, the self-compaction is supposed to be much more efficient and intensive than elsewhere due to the much higher snow/ice level at the end of the northern fall (compare Figure 5 to Figure 4). Once the low-density fresh snow compacts into denser porous slab ice, sintering of the ice crystals can happen at an accelerated pace due to positive feedback (Eluszkiewicz et al., 2005). Furthermore, much smaller amounts of snowfall during the early to-middle . Individual stretching is applied to each sub-frame to emphasize local contrast. Refer to Figure 5 for the corresponding level variation time series from Mars Orbiter Laser Altimeter laser altimetry. The accentuated texture at L s 1° points to the existence of translucent slab ice, which can be due to snow densification during the early-to-middle northern winter as indicated by the off-season decrease of nearly 2 m. Then, dark dust fans and polygonal cracks begin to form along the crestlines and sun-facing slip-faces due to cold-jetting mechanism (L s 35° and L s 62°). Translucency of the slab ice enables the sunlight to penetrate through and sublimate the dry ice at the substrate. Eventually, pressure builds up, cracking the slab ice and erupting the dust materials to the surface which are then blown away by the winds to form the dust fans. The sublimation of the Seasonal North Polar Cap happens quickly during the middle-to-end of the northern spring. Thin layer of water ice persists to the beginning of the northern summer (L s 91°), some bright materials between the dunes, most likely gypsum and other evaporitic minerals, can linger throughout the northern summer (L s 114° and L s 174°).
northern winter compared to that in the middle-to-late northern fall and late northern winter were observed at the northern hemisphere, especially at latitudes higher than 80°N (refer to Figure 6 in Gary-Bicas et al. (2020) that shows the CO 2 cloud opacity which is most likely associated with polar snowfall). This decrease in snowfall can be related to the "solsticial pause" in low-altitude transient planetary wave activity (Lewis et al., 2016), which can inhibit poleward transport of dust and water ice particles that act as condensation nuclei for atmospheric CO 2 . The phased accumulation can then accentuate the metamorphism effects during the early-to-middle northern winter. The metamorphism combined with a self-cleaning process, in which dust grains burrow down to the ground or released through the upper surface (Kieffer, 2007), can ultimately result in translucent slab ice as can be observed at L s 1° in Figure 12. It should be noted that the earliest HiRISE images we examine within Olympia Undae (around L s 344° at the end of northern winter) have already featured translucent slab ice with low surface reflectance. Unfortunately, there are no optical images available in the early to middle northern winter to confirm the existence of translucent ice due to the dark polar night. Even we have access to those images, the seasonal cap can be obscured by the polar hood of optically thick water ice clouds (Benson et al., 2011;Brown et al., 2016;Calvin et al., 2015;Navarro et al., 2014). The limited accumulation in the late northern winter as can be observed in Figure 5 could also be the result of this phased snowfall at high latitudes (Gary-Bicas et al., 2020).
The off-season increases during northern spring mainly happen earlier than L s 30°. Unfortunately, no MOLA profiles were acquired during this period due to the solar conjunction. However, the assumption of large quantities of snowfall during polar cap recession at springtime within Olympia Undae would be unlikely, as both MOLA (MY24-25) and MCS (MY29-33) measurements have indicated a declining cloud opacity starting from L s 330° which drops to zero at L s 360° (the end of northern winter; Ivanov & Muhleman, 2001;Gary-Bicas et al., 2020). Meanwhile, tracking of the seasonal polar cap edge using thermal infrared and spectroscopic measurements has indicated a nearly monotonic retreat during the springtime and CO 2 snow/ice is almost completely absent from the surface of Olympia Undae by L s 30° (Appéré et al., 2011;Kieffer & Titus, 2001;Piqueux et al., 2015).
The translucence of the slab ice observable in HiRISE images in late northern winter and early northern spring is consistent with the phased accumulation observed in high-latitudes at the north pole. Relatively low snowfall during early-to-mid northern winter combined with gravity compaction, re-crystallization, and dust burrowing down the seasonal layer would mean the seasonal deposits can become increasingly translucent during this period. Increased snowfall and formation of a water ice veneer thereafter can gradually lower the translucence ( Figure 12). In light of the evolution of translucency of the seasonal deposits, an additional tentative explanation for the off-season anomalies within Olympia Undae is that the increasing translucency of the ice/snow deposits during early-to-mid northern winter enables the MOLA laser pulses at 1,064 nm to partially penetrate, biasing the laser one-way ranges and lowering the level measurements. Then, as the translucency decreases during later northern winter and early-to-mid northern spring, MOLA gradually re-measures to the surface, recovering the level estimates. However, it is concluded that while the translucent slab ice can let in visible sunlight, it remains largely opaque in the thermal infrared spectrum (e.g., Chinnery et al., 2020). For that, Chinnery et al. (2020) acquired a mean penetration depth of 4.7 cm into the slab ice for wavelengths between 300 and 1,100 nm. However, this experiment was conducted in a state (grain size/cracks) that may be not representative of Martian ice. We also carry out our own modeling using a detailed and comprehensive radiative transfer model for rough CO 2 slab ice, as described in Andrieu et al. (2015). Our preliminary results demonstrate that specular reflection at the atmosphere-ice interface does not dominate the return signal. A majority of the pulse energy (60%) can manage to penetrate the slab ice. However, how these would introduce biases to the MOLA heights depends on the interaction of the returned waveforms with low-filtering channels and prefixed pulse energy thresholds of MOLA for measuring the time-of-flights. These are complicated to quantify and require extensive simulations on how the laser pulses interact with the sloped slab ice within Olympia Undae. Meanwhile, detailed analysis on how the simulated returned laser waveforms interact with MOLA measuring characteristics needs to be performed in the future. If the penetration theory is indeed behind those unexpected features observed at Olympia Undae, then this can also have implications for other laser altimeter missions, including those to the icy moons of Jupiter and Saturn, for example, the Ganymede Laser Altimeter (GALA; Hussmann et al., 2019).

Asynchronous Extent and Volume Accumulation
The volumes of both seasonal polar caps peak at approximately the end of winter (refer to Figure 3). Consistently, they continue to gain mass until late winter (e.g., Karatekin et al., 2006;Kelly et al., 2006;Schmidt & Portyankina, 2018). These are in contrast to the fact that the aerial extents of both seasonal polar caps maximize at the very beginning of winter (around L s 5° after the winter solstice) and start to shrink thereafter (Piqueux et al., 2015). The asynchrony can be explained by the fact that while CO 2 snow/ice at the lower-latitude polar cap edge begins to sublimate earlier in the wintertime, CO 2 condensation/deposition is still happening at a faster rate toward the higher-latitude colder polar regions. This process can continue throughout the dark winter polar night.

Bulk Density Evolution of the SNPC
To examine the temporal evolution of the bulk density of the SNPC as we have done for the SSPC (Xiao, Stark, Schmidt, et al., 2022), mass change of the SNPC measured from gravity variations (Karatekin et al., 2006) and its corresponding volume during MY24 are compared ( Figure 3). Volume time series of the SNPC is fitted using the Gaussian Process regression. Then, volume measurements are sampled from the fitted curve with an interval of 2° in solar longitude from L s 220° in MY24 to L s 450° in MY25. The bulk density of the SNPC is then obtained from dividing its mass by the corresponding volume measurement (black line and crosses in Figure 13). Uncertainties of the derived bulk density are also obtained by evaluating the lower and upper bounds of the mass measurements from Karatekin et al. (2006) against that of the volume measurements treated as the 95% confidence bounds from the Gaussian Process fitting. The general evolution of the bulk density is in contrast to the normal assumption that the bulk density should continuously increase throughout fall, winter, and spring (e.g., Matsuo & Heki, 2009), as also illustrated in Xiao, Stark, Schmidt, et al. (2022) at the south pole. Thus, nonlinear modeling of the bulk density needs to be incorporated into future Mars volatile circulation models. Despite continuous increase of the mass from Karatekin et al. (2006) throughout the northern fall and winter, we can see phased increases of the volume as in middle-to-late northern fall and middle-to-late northern winter, respectively. These phased volume accumulation, and resultant unconnected phases of decreasing bulk density, can be due to the phased snowfall, especially at high latitudes (Gary-Bicas et al., 2020). While between these two volume increase phases, the amount of snowfall is relatively small, so metamorphism can play a main role which increases the bulk density. Stabilizing volume during most of the northern spring in MY25 is considered as an artifact (yellowish shade in Figure 13) as the SNPC is completely sublimated away by the end of northern spring, as indicated by its area measured by thermal spectrometers (Piqueux et al., 2015) and its mass by gravity variations (Karatekin et al., 2006). As a result, the decreasing trend of the bulk density after a transient increase during the early northern spring is most likely invalid.
We should note that, unlike for the SSPC, the existing mass estimates of the SNPC are largely inconsistent with maximums ranging from 3.5 to 6 × 10 15 kg. Thus, the temporal bulk density analysis for the SNPC here can only stand if the mass estimate from Karatekin et al. (2006) is correct. Indeed, considering a lower bound of 3.5 × 10 15 kg and a upper bound of 6 × 10 15 kg, the bulk density for the SNPC when its volume maximizes at 9.6 × 10 12 m 3 can be calculated at 365 kg/m 3 and 625 kg/m 3 , respectively. The lower bound of 365 kg/m 3 is significantly lower than the corresponding value of 660 kg/m 3 at the south pole (Xiao, Stark, Schmidt, et al., 2022) and even slightly lower than the threshold of 420 kg/m 3 for snow (≥74% porosity; Mount & Titus, 2015). This comparison can indicate higher quantities of atmospheric depositions at the north pole than at the south pole, which is consistent with the surface emissivity measurements from MCS measurements (Gary-Bicas et al., 2020) and grain size index from Thermal Emission Spectrometer (TES; Haberle et al., 2004). Meanwhile, the upper bound is just slightly lower than that of the SSPC.
Unfortunately, the mass change rates used in Xiao, Stark, Schmidt, et al. (2022) are not available for the north pole to constrain the bulk density of the snow/ice deposits at local scales. The non-negligible presence of water ice and dust within the seasonal polar cap in the north compared to that in the south can significantly modify the  Figure 3, and the 95% confidence bounds (blue circles and shade). Also shown are the mass evolution of the SNPC from time-variable gravity measurements (red stars and shade; Karatekin et al., 2006), and the resultant evolution of its bulk density (black crosses and shade). For clarity, only every four markers are shown (a temporal gap of 8° in solar longitude). Yellowish shade denotes the period in northern spring during which the measured volume of the SNPC and its resultant bulk density are invalid (Section 5.2). The bulk density threshold of CO 2 snow (≤420 kg/m 3 ; Mount & Titus, 2015) is marked for reference. surface thermal emissivity and albedo of the seasonal snow/ice deposits, which need to be accurately modeled in the energy flux in order to obtain reliable mass change measurements (Schmidt et al., 2009(Schmidt et al., , 2010. Thus, it is not feasible to spatially resolve the bulk density of the SNPC as has been done to SSPC (Xiao, Stark, Schmidt, et al., 2022).

Conclusion
We retrieve seasonal level variations of Martian SNPC by applying a local co-registration strategy complemented with a post-correction procedure based on pseudo cross-overs to MOLA profiles. These techniques have been proposed in Xiao, Stark, Steinbrügge, et al. (2022) and previously been applied to the SNPC at the south pole (Xiao, Stark, Schmidt, et al., 2022). The level variation time series is generated in individual grid element of 0.5° in latitude from 60°N to 87°N and 10° in longitude. The maximum level over the Residual North Polar Cap is just 1.3 m as compared to 2.5 m at the south pole. Exceptionally, maximum level in the dune fields at Olympia Undae can be up to 3.8 m. Off-season decreases up to 3 m during the northern winter at Olympia Undae are observed, which could possibly be attributed to snow/ice metamorphism, accentuated by reduced snowfall during this period. At the same time, off-season increases of up to 2 m during the northern spring are observed. Theoretically, these off-season anomalies may also be attributed to artifacts in MOLA heights due to interaction of the laser pulses with the dynamic properties (footprint-scale slope, roughness, and albedo) of the dune surface, or related to the penetration of the laser pulses into the translucent slab ice. We plan to examine these hypotheses in dedicated future studies. Early height accumulation during beginning winter is much higher and more concentrated as a whole at the north pole than at the south pole. The peak volume of the SNPC at the end of northern winter is calculated at around 9.6 × 10 12 m 3 , which approximately equals that of the SSPC (9.4 × 10 12 m 3 ). The average bulk density of the SNPC when its volume maximizes is constrained to be in the range of 365 to 625 kg/m 3 . In addition, the bulk density of the SNPC can go through phased decreases due to phased snowfall at northern high-latitudes.
This study examines both the temporal and spatial level evolution of the SNPC, which should be incorporated as constraints in future Martian climate and volatile models (Becerra et al., 2021). As the next step, we will try the SHAllow RADar (SHARAD) radar altimetry to cross-validate the MOLA measurements, and to derive the long-term seasonal level evolution of the seasonal polar caps of Mars (e.g., Raguso & Nunes, 2021;Steinbrügge et al., 2021). These results can be adopted to verify the observed off-season decreases and increases at Olympia Undae, and to examine if they can annually repeat themselves in the long term. Besides, the long-term records will also be important for assessing the stability of the Martian residual polar caps which is still debated. Proposed laser altimeters dedicated to Mars' polar mapping, such as those instruments included in the Orbital Radio science and Altimetry for CLimate Experiment mission (Genova, 2020), the Atmospheric/Surface Polarization Experiment at Nighttime lidar (Brown et al., 2015), and a proposed next-generation Mars polar explorer (Oberst et al., 2022), are expected to gain further insights into seasonal evolution of snow/ice levels at Martian poles.

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