Mercury’s Hidden Past: Revealing a Volatile-dominated Layer through Glacier-like Features and Chaotic Terrains

The discovery of global elemental volatile compositions, sublimation hollows, and chaotic terrains has significantly reshaped our understanding of Mercury’s geology. These findings suggest the existence of volatile-rich layers (VRLs) extending several kilometers in depth, challenging the traditionally held view of a predominantly volatile-devoid Mercury crust. However, the precise nature and origin of these VRLs remain to be elucidated. The Raditladi basin exhibits morphologies analogous to terrestrial and Martian glaciers. These geomorphological features are potentially derived from impact-exposed VRLs, likely constituted of halite, other semivolatile salts, or organic volatiles. The distinctive rheological traits of substances such as halite substantiate this hypothesis. The inference posits a potential ubiquity of VRLs on a planetary scale, albeit potentially ensconced at considerable depth in specific regions. North polar chaotic terrains elucidate the VRLs’ genesis and temporal evolution. The intense fragmentation of heavily cratered landscapes during their formation indicates a composition dominated by volatiles. This finding postulates a phase of volatile-enriched crustal accretion predating the Late Heavy Bombardment (∼3.9 Ga). Regardless of lost mass, the unaltered basal elevation post-collapse signals a transition to a volatile-free stratum. The exposure of an exhumed lithological substrate within Mercury’s stratigraphy, identifiable in gravimetry as an impacted paleosurface, contests the magma ocean differentiation concept for VRL formation. It infers a grand-scale construct originating from depositional processes, possibly due to the collapse of a transient, hot primordial atmosphere.

While these elements were detected in a superficial layer roughly between 10 cm and 100 μm thick (Lawrence et al. 2013;Neumann et al. 2013;Deutsch et al. 2017;Chabot et al. 2018), the finding of potential sublimation pits, also known as hollows (Figure A1), sheds light on the distribution of materials within Mercury's volatile-rich subsurface.These hollows, identified as flat-floored, rimless depressions featuring bright interiors and halos (Blewett et al. 2011(Blewett et al. , 2013(Blewett et al. , 2016;;Helbert et al. 2013), have an average depth of 24 ± 16 m, interpreted by Blewett et al. (2016) to indicate that volatiles exist as a pervasive constituent extending well beneath the planet's surface.
The Caloris Basin Antipodal Chaotic Terrain (CBACT) and several other recently analyzed chaotic terrains manifest significant elevation drops ranging from a few hundred meters to several kilometers (Rodriguez et al. 2020).Such observations imply that the constituents of the ancient, cratered landscapes have experienced disaggregation, possibly triggered by the sublimation of substantial volatile-rich layers (VRLs; Rodriguez et al. 2020).Considering a possible inception of the CBACT ∼3.9 Ga and its primary developmental cessation ∼1.8 Ga, it is plausible that this disaggregation process spanned a significant duration (Rodriguez et al. 2020).Wright et al. (2020) interpreted the degradation patterns within the Caloris ejecta blanket as indicators of devolatilization over large areas of impact-overturned VRL materials, also underscoring the theory of a stratigraphy composed of thick VRLs.
In the present study, we propose that VRLs exhibit a broader spatial distribution than previously understood, encompassing significant subterranean regions.Our findings indicate that these concealed areas are revealed in impact craters, unveiling the intriguingly widespread occurrence of VRLs across the planetary surface.These excavated materials gave rise to morphologic characteristics reminiscent of glaciers and especially glacier-like Martian lobate debris aprons.Moreover, through our meticulous examination of polar chaotic terrains, we assert that the VRLs embody substantial condensation features, potentially formed prior to the Late Heavy Bombardment (LHB) due to an early atmospheric collapse.We propose that these materials experienced planet-wide accumulation and regionally extensive burial.

Geomorphologic Mapping Approach
Our investigation focuses on analyzing two regions: (1) the Raditladi impact basin, which has a diameter of ∼263 km and is situated at 27°N and 241°W (Figure 1).According to Strom et al. (2008) and Prockter et al. (2009), the basin likely formed ∼1 Ga; and (2) a chaotic terrain located in the northern polar part of the planet, Borealis Chaos (BC), which is roughly centered at 76°N and 150°E (Figure 1).
We used ArcGIS software from the Environmental Systems Research Institute (ESRI 2016) to compile and geo-register the relevant MESSENGER MDIS data products.The geo-registration process was carried out manually in ArcGIS, whereby we tied the locations of features to the base map and rectified the image through rubber-sheeting based on the tie-point locations.Detailed information about ESRI can be found at http://www.esri.com.
We conducted geological feature mapping, which included the lobate aprons, hollows of the Raditladi impact basin, and various polar chaotic terrain types.To achieve this, we used two sources of data: the MESSENGER MDIS global mosaic base map with a resolution of 166 m pixel −1 and the MESSENGER global Digital Terrain Model (DTM) version 2 from 2016 October, with a resolution of 665 m pixel −1 and ∼1 m vertical precision.
In addition, we used the MESSENGER Regional Targeted Mosaics (RTMs) with a resolution of 15-24 m pixel −1 , which we obtained and geo-registered in the Raditladi basin.Moreover, we employed in-house DTMs derived from stereophotoclinometry (SPC) with 20 m ground sample distance (GSD) resolution to perform detailed morphological and topographical characterizations of the case study areas (Figures A7-A9; see details on SPC methodology in Appendix A).

Raditladi Impact Basin
Upon detailed examination of the Raditladi Basin, distinct hollow clusters are identified, specifically situated on its peak ring and the adjacent lobate depositional aprons (Figure 2(A)).These aprons cover an extensive expanse, nearly 300 km of the overall 450 km circumference of the peak ring, and they clearly are associated with the peak ring mountains (Figure 2(A)).Further, their cross-sectional morphology exhibits convex upward forms, terminating abruptly at steep ends (Figures 2(B), (C); elevation profiles X-X′ and Y-Y′, respectively).
The lobate formations consistently extend to runout distances of tens of kilometers from the peak ring scarps, in some instances superimposing extensional faults on the adjacent floor (Figures 2(A), (C)).Several hollows exhibit depths representing a significant portion of the overall lobes' thickness.As an example, in Figure 2, profile Y-Y′ showcases a hollow approximately 75 m deep, a depth three times the average of 24 ± 16 m, as specified by Blewett et al. (2016), nested within a deposit locally measuring around 250 m in thickness.Further, extensive hollow clusters are discernible both on and along the peak ring, with a portion merging to form relatively flat areas on the peak ring (Figure 3(A)).
Numerous moated craters manifest across multiple lobes, hallmarked by rim areas exhibiting depressions that shape semicircular troughs (Figures 2(A), 4(A)).Contrasting with the peak ring lobes, the basin's interior rim flanks project depositional aprons resembling rotational landslides and are devoid of superimposed hollows (Figure 5).
Various lobes exhibit a significant presence of moated craters, marked by semicircular troughs along their rim areas (Figures 2(A), 4(A)).Contrasting the peak ring lobes, the interior rim flanks of the basin feature depositional aprons reminiscent of rotational landslides, notably without any superposed hollows (Figure 5).Another impact peak ring has been identified within the 125 km diameter Eminescu crater (10°.66N;245°.79W; Figure 6).This peak ring is surrounded by widespread depressions that present outlines mirroring those of lobes in the neighboring Raditladi's peak ring.The Eminescu crater, akin to the Raditladi basin, is estimated to have formed ∼1 Ga (Schon et al. 2010).

Borealis Chaos
Borealis Chaos (BC; delineated in Figures 1 and 7) represents a region of significant geological complexity, predominantly localized within the boundaries of Borealis Planitia, Mercury's expansive northern plains.This area extends longitudinally from approximately 90°E to 240°E, while latitudinally, it reaches from the north pole to approximately 50°N (Figure 1).Notably, the portion of the region situated longitudinally between 150°E and 240°E encompasses previously detected high gravity anomalies (Qingyun et al. 2018) (Figure 7).
BC is characterized by the landscape segmentation of the regional heavily cratered terrains.This segmentation is such that craters across a range of diameters exhibit consistent degradation and discontinuity in their rims, as denoted by arrows 1-3 in Figure 8(A) (spanning diameters approximately from 10-160 km).The planar component of BC displays a nearly uniform topographic base, established at an elevation of roughly −4000 m, as exhibited in the elevation profile in Figure 8.A close-up view reveals a crater-saturated surface with some highly degraded craters that have merged into irregular depressions (e.g., feature 1 in Figure 8(B)).The terrain also includes larger irregular depressions (e.g., feature 2 in Figure 8(B)) and broad plains harboring mesas (e.g., feature 3 in Figure 8(B)).

Volatile-rich Mass Flows and Landscape Shaping in the Raditladi Basin
Our detailed examination of the Raditladi Basin unveils distinct distributions of hollows (Figure 2(A)), suggesting a preponderance of volatile-rich materials within the peak ring that were likely derived from uplifted shallow crustal zones (Kring et al. 2016).The substantial depth of hollows overlaying the lobes-locally reaching up to a third of their total thickness -hints at a composition primarily characterized by volatiles (profile Y-Y′ in Figure 2).This observation, coupled with the deposit's morphology-marked by lobated fronts and a convex upward surface topography with abrupt termini (Figures 2(B), (C), elevation profiles X-X′ and Y-Y′)-mirrors the typical behavior of viscous material, thereby prompting us to propose a glacier-like origin.
Furthermore, the seemingly flat ridge tops within the peak ring materials (Figure 3(A)) suggest a demarcation differentiating volatile-rich compounds from nonvolatile rocks.Our hypothesis proposes these lobes could have originated from the flow of exposed VRLs, a process that likely contributed to reducing the peak ring's relief A fundamental precondition for the preservation of the lobes involves the development of a refractory lag capable of arresting or significantly mitigating devolatilization.The manifestation of such a lag corresponds with the observed prevalence of moated craters (Figure 4(A)), distinguished by their rim regions yielding semicircular troughs.We propose that these moats were formed due to the impact-induced uplift, exposure, and consequent erosion of the volatile-rich substrates sequestered beneath the refractory lag along the crater rim zones (Figure 4(B)).
Following this hypothesis, the initiation of hollows might correspond with zones where the refractory layer-originating  N01_003818_2236022). Detailed procedures for RTM processing can be found in the supplementary materials, with context and geographic location provided in Figure 2(A).(B) Cross-sectional reconstruction, drawn along line M to M′, postulates the formation of moats due to impact-induced exposure of near-surface, volatile-rich materials along the edges of some small (subkilometer to a few kilometers) craters.(C) The modeling of surface solar heat waves reveals a relatively shallow penetration depth (∼3-8 m), consistent with the long-term retention of volatiles (see Sections 4.2.2.4 and 4.2.2.5).The presented profiles represent the thermal extremes, with the warmest portion of daytime (upper curves) and the coldest point of nighttime (lower curves).This finding, aligning with the hypothesis that hollows are recent geological features, extends the same process to explain the formation of moated craters.The model incorporates a range of parameters, including a peak solar irradiance of 9000 W m −2 , emissivity set at 0.95, an albedo of 0.06, and a geothermal flux of 20 mW m −2 .Thermal conductivities are specified as K t = 3 W m −1 K −1 for intact rock, K t = 1.75 W m −1 K −1 for a composite substance, and K t = 0.3 W m −1 K −1 for highly porous soil.from the desiccation of the surfaced VRL-was disturbed or fractured owing to impact events or other geological dynamics such as faulting.
Expanding on our hypothesis, we posit that the discernible lobate depressions adjacent to the Eminescu crater's peak ring (Figure 6) were likely formed due to the near-total sublimation and eradication of glacier-like deposits, potentially enabled by a comparatively thin refractory layer.In addition, these observations underscore the metastable condition of Mercury's glacier-like features, thereby suggesting that older impact craters may have once hosted such features that have since vanished.Consequently, their absence in older craters does not inherently indicate a discontinuity in VRLs at depth.
It is pivotal to underscore that these deposits are interpreted solely as morphological analogs.Considering Mercury's exospheric properties, these glacier-like features are not hypothesized to have resulted from atmospheric mechanisms and are not considered to have involved H 2 O ice but are perceived as outcomes of endogenic processes acting on materials that are volatile at Mercury surface conditions.

Lobate Debris Apron Analog Considerations
The lobes observed in the Raditladi Basin manifest striking similarities to Martian Lobate Debris Aprons (LDAs), interpreted as ice glaciers covered by debris on Mars (Mangold 2003;Arfstrom & Hartmann 2005;Hauber et al. 2008;Baker et al. 2010;Séjourné et al. 2019).This comparison lends credence given the close match in surface gravitational acceleration between Mars and Mercury-3.71m s −2 and 3.7 m s −2 , respectively-a key factor in determining typical driving stresses for LDA propagation.
The Raditladi lobes and Martian LDAs (Figures A2(A), (B)) both exhibit restricted runout distances (barely spanning tens of kilometers from their originating scarps; Figures 2(A Squyres 1978;Mangold & Allemand 2001;Pierce & Crown 2003;Chuang & Crown 2005, 2009;Li et al. 2005;Berman et al. 2015;Schmidt et al. 2019).When feasible to quantify, these reveal similar maximal thicknesses-typically fluctuating between approximately 200 and 250 m.These dimensions, apparent in Figures 2(C) and (D) for the Mercurian lobes under study, align with the observations made by Schmidt et al. (2019) for Martian LDAs.The consistency in these characteristics underscores our supposition that the flow dynamics governing these geomorphic feature emplacements on Mars and Mercury share fundamental resemblances.
The compelling morphological parallels between the Raditladi features and Martian LDAs suggest shared volatile-rich origins.Our analysis posits that these glacier-like deposits in Raditladi likely owe their existence to a volatile-enriched subsurface stratum that became exposed in some sloping areas.Similarly, Martian LDAs are hypothesized to arise from a vast, ice-dense layer embedded within the Martian highlands (Pierce & Crown 2003;Kargel 2004;Chuang & Crown 2009).
The Raditladi lobes and Martian LDAs demonstrate significant parallels in relation to the preservation, flow, and partial loss of volatiles.Notably, the presence of hollows throughout the Raditladi lobes (Figure 2) and the evidence of relatively recent sublimation pits on Martian LDAs' surfaces A close-up of BC's floor reveals a crater-saturated surface, including notably degraded craters morphing into irregular depressions (feature 1).Larger irregular depressions (feature 2) and even more extensive plains with enclosed mesas consisting of heavily cratered plateau residues (feature 3) can also be seen.A regionally consistent topographic base around −4000 m indicates the interface between a regional VRL and the gravity-detected cratered, lava-covered paleosurface (refer to Figure 7).Data source: both panels consist of a MESSENGER Mercury global DEM (∼665 m px −1 ) draped over a MESSENGER MDIS global base map (∼166 m px −1 ; courtesy of NASA/ JPL/GSFC).(Pierce & Crown 2003) indicate a long-term subsurface sequestration and conservation of volatiles, bolstered by observations of locally exposed, radar-identified debris-covered ice on Mars (Plaut et al. 2009).
Furthermore, the stratigraphic similarities between the Raditladi lobes and Martian LDAs go beyond shared volatile origins and retention mechanisms, encapsulating the properties of their corresponding nonvolatile topmost stratigraphic components.The presence of moated craters within the lobes (Figures 2(A), 4(A)) attests to a robust lithic layer that encapsulates and stabilizes underlying volatiles (Figures 4(B), (C)), an attribute mirrored in Martian geophysical data (Kress & Head 2008).
Further, the origins of the lobes on Mercury via lowviscosity halite flows find a terrestrial parallel in Iran's diapirsourced salt glaciers (Talbot & Pohjola 2009).There are critical morphological and geological similarities between Mercury's lobes and these terrestrial salt glaciers, which include: On Mars, the genesis of LDAs, commonly linked to the expulsion of substances from geological formations enriched with ice, characteristically spread across broad plains encircling their origin area (Pierce & Crown 2003;Kargel 2004;Chuang & Crown 2009).On the other hand, terrestrial glaciers primarily originate from the gravitational compaction and downslope movement of accumulated and recrystallized snow within geographically restricted valleys, resulting in a distinctive geomorphological footprint defined by channeled glacial paths and the descent of valleys.
However, a terrestrial geomorphological feature that offers a persuasive counterpart to Martian LDAs resides in the salt glaciers of the Zagros Mountain range (Figure A2(C)).These captivating landforms, which emanate from salt diapirsstructures tectonically akin to the cryospherically enriched block formations seen on Mars-showcase parallel morphodynamic characteristics.This makes them a notably robust comparative framework for Martian LDAs.

Comparative Analysis of Volatile Element Abundance in Mercury, Mars, and Earth
Another intriguing parallel between Mercury and Mars arises in the context of their volatile constituents.Specifically, Mercury's moderately volatile elements-namely Na, S, Cl, and K-have been found to occur in proportional abundances that echo those observed for Mars (Nittler et al. 2018).In addition, an inspection of the ratios of potassium-to-thorium (K/Th) and chlorine-to-potassium (Cl/K), as derived from MESSENGER's Gamma-Ray and Neutron Spectrometer (GRNS) readings, reveals that these quantities markedly exceed terrestrial abundances, drawing a closer alignment with the elemental constitution of Mars (Peplowski et al. 2012;Nittler & Weider 2019).
This divergence from terrestrial norms becomes even more pronounced when focusing on specific elements.Cl concentrations in bulk-Earth crust plus mantle are around ∼36 ppm (Kargel & Lewis 1993).Yet, on Mercury, these concentrations dramatically rise, reaching ∼1200 ppm near the equator and soaring to ∼2500 ppm in the high northern latitudes (Evans et al. 2015;Nittler & Weider 2019).A similar trend emerges for S. While terrestrial mantle plus crust contains a mere few hundred ppm of S, Mercury's surface exhibits an average abundance of ∼4 wt% (Nittler et al. 2018).

Compositional and Viscosity Constraints on the Glacierlike Features through the Simulation of Impact and Flow Processes on Mercury
The discernment of potential glaciogenic flow traits within our observational scope suggests the manifestation of volatile constituents with capabilities for transient ductility.Subsequently, we quantitatively assess the possible compositional makeup and viscosity associated with the flow material prevalent during the formation of Mercury's glacier-like deposits.

Modeling the Raditladi Basin Post-impact Stratigraphy and
Cooling History

Impact Excavation Modeling Approach
The development of the Raditladi crater, a peak-ring geological formation with a terminal rim-to-rim diameter (D r ) of approximately 250 km, was simulated employing the iSALE2D hydrocode (https://isale-code.github.io/).The transient crater diameter was estimated using the methodology developed by Abramov & Kring (2005): In the equation above, D r represents the final rim-to-rim crater diameter, D tr signifies the transient rim-to-rim crater diameter, and D Q denotes the transition diameter from simple to complex craters (which stands at 10.3 km m for Mercury, according to the gravitational scaling by Croft 1985).By utilizing this relationship, a D r of 250 km equates to a transient rim-to-rim crater diameter (D tr ) of 140 km and a transient crater diameter at the pre-impact surface (D tc ) of 115 km (where the corresponding radii are related by R tc = R tr /1.2; based on Melosh 1989).This last value, D tc , is then used in conjunction with the pi-group scaling laws (Schmidt & Housen 1987) to deduce the projectile diameter, assuming a standard impact velocity of 30 km s −1 for Mercury.
The development of Raditladi is simulated as a vertical collision by a dunite sphere with a derived 10 km diameter.The model operates at a spatial resolution of 400 m, with the target surface portrayed by a 435 × 380 element grid.The initial temperature is calibrated to 440 K (163 °C), with no geothermal gradient incorporated.In the primary model, the structural profile (Childs 1993), from bottom to top, is composed of a 40 km dunite mantle, a 30 km basalt crust, and a 2.8 km VRL characterized by calcite (Figure 9(A)).These materials are representative of an envisioned structure, but do not explicitly demonstrate compositional constituents.
Additional models incorporate variations such as a dunite crust, an extended 11 km VRL, and/or a regolith overburden of 1.2 km situated over the VRL (additional details are provided in supplementary materials).

Modeling Thermal History Approach
Utilizing the versatile, finite-difference heat transfer program HEATING 7.3, developed by the Oak Ridge National Laboratory, we simulated the post-impact cooling phase of the Raditladi crater.The preliminary temperature distribution was obtained from the final stage of the iSALE simulation that incorporated a basaltic crust.This model employed thermal and physical parameters specific to basalt, as provided in the HEATING materials library (Childs 1993).Our model, which presents a 2D radial symmetry (HEATING geometry type 3, rz cylindrical), used a 435 × 380 element grid with a spatial resolution of 400 m per element while integrating radiative cooling at the surface boundary.
In the simulated model, a VRL represents the uppermost stratigraphy (Figure 9(A)), a setup that aligns with mapping data demonstrating a lack of hollows in the basin's central region and their sporadic distribution in the plains interspersed between the peak ring and crater rim (Figure 2(A)).Our model, Figure 9.This illustration presents a perspective view of the Raditladi impact crater basin.The left half of the image reveals the subsurface temperature distribution right after the impact, as projected by the iSALE2D hydrocode simulation of the crater's formation.In contrast, the right half reveals the stratigraphic makeup of the impact-excavated upper section, comprising a 2 km thick VRL surmounting a multi-kilometer-thick volcanic crust.Note that the peak ring region (highlighted in red boxes in panels (A) and (B)) emerges as a distinctive area, marked by a synergy of high relief and the outermost confines of the impact-triggered thermal disturbance (as further detailed in simulation close-ups (B) and (C)).
though simplified and limited in reflecting complex processes of peak ring formation, is consistent with the origin related to the uplift of regionally buried VRL materials (Kring et al. 2016).

Impact Excavation Modeling Results
The output from the iSALE2D hydrocode model chronologically renders a sequence of post-impact alterations (Figure A3).The maximum impact depth is achieved at the 30 s post-impact interval.Inverted stratigraphy within the ejecta curtain becomes discernible at the 170 s post-impact interval.At 300 s post-impact, the initial emergence of a transient central peak, along with the formation of an ejecta blanket, is detected.The basin stabilizes into its final form approximately 700 s after the impact event (Figure 9(A)).
The modeled stratigraphy implies the potential removal of a regional VRL within the peak ring area due to the impact, revealing the basal basaltic layer (Figure A3).It also suggests a possible, albeit diminished, preservation of the VRL layer within the annular zone bounded by the peak ring and the crater rim (Figure 9(A)).For potential alternative stratigraphic arrangements, see Figures A4-A6.Each configuration suggests retention of the VRL at the peak ring site.
In an in-depth analysis of the post-impact thermal dynamics, we implemented a comprehensive simulation to chart the thermal trajectory of the basin over a timescale of 1 Myr following the impact.The results indicate that the substantial thermal energy resulting from the initial impact was a key factor in increasing the temperature of the peak ring, situating it within a thermal band approximately spanning 300°-500 °C (Figures 9(B) and (C)).

Computational Modeling of Mercurian Glaciogenic Flow Features
Subsequently, we present our computational flow models meticulously designed to decode deposition mechanisms associated with Mercury's postulated glaciogenic structures.These models are precision-tuned considering the thermal architecture inferred from previous sections and an assessment of feasible volatile constituents.

Halite: A Probable Constituent
The exact composition of the glacier-like flow materials, however, remains undetermined.In the impact model, we modeled the VRL as calcite.Here, we consider another substance, halite, for the flow material, which would be a component of the VRLs.Both minerals are abundant and widespread in the solar system and are not mutually exclusive.Preliminary examinations of hollows suggested the regional existence of various substances near Mercury's surface, including chalcogenides, pnictides, sulfosalts, and different salts (Kargel 2013).Furthermore, Rodriguez et al. (2020) singled out halite and other mildly volatile salts as potential constituents, guided by their comprehensive assessment of cosmochemical elemental abundances.
These considerations collectively substantiate the plausible occurrence of slightly volatile minerals on Mercury, specifically halides-most notably sodium chloride (NaCl, or halite) -composed of elements such as sodium (Na), potassium (K), chlorine (Cl), and sulfur (S).These minerals may not merely be present but could potentially be the dominant constituents within the region.Rodriguez et al. (2020), building on Kargel (2013), underscored that halite, owing to its somewhat volatile characteristics, might incite terrain degradation over geological timescales, thereby possibly facilitating the genesis of chaotic terrains.
It is noteworthy that halite and other halide salts are recognized as substantial near-surface components on Earth and Mars.Adding weight to this, recent confirmations of halite and hydrated NaCl on Ceres (De Sanctis et al. 2020;Bramble & Hand 2022), in addition to their presence within certain meteorite types relevant to Mercury (Rodriguez et al. 2020), make the postulation of halides existing on Mercury even more credible.Notably, halite shifts to blue coloration due to radiolytic alterations (Zelek et al. 2014).This phenomenon might potentially explain the observed comparatively bright blue hue of Mercury's hollows (Munaretto et al. 2023).
Indeed, sodium compounds display significant volatility, especially when combined with sulfur or chlorine, forming Na 2 S or NaCl (halite).These compounds could potentially undergo extensive sublimation on Mercury, similar to that observed with water ice on Callisto.On Callisto, sublimation has mediated considerable modifications to cratered terrains over geological timescales spanning billions of years (see Figure S14 in Rodriguez et al. 2020).
In our simulations, we propose sulfur (S 8 ) and NaCl as the likely predominant constituents of the glacier-like features.
Here, S 8 denotes the most stable allotrope of sulfur, which consists of a cyclic arrangement of eight sulfur atoms.Based on the initial temperature conditions within the peak ring, as derived from the impact simulation (Figure 9), it is evident that the temperature ranges are insufficient to induce the melting of sulfides such as MgS, CaS, and CS 2 .However, these ranges intersect with the melting point spectrum of sulfur (S 8 ), specifically between 110 °C and 440 °C (Figure B1).
Therefore, sulfur, which could have originated from the breakdown of sulfides, represents a probable volatile linked with the formation of hollows (Phillips et al. 2021) and hence presents itself as a viable candidate (Figure 3(E)).While other volatiles, such as organic compounds, specifically stearic acid, and fullerenes, have also been suggested (Phillips et al. 2021), they are not considered in this study due to the lack of compelling evidence supporting their substantial generation and presence.
Our simulations advance our grasp of halite's rheological properties, demonstrating ductility in contrast to S 8 (Figure 3(D)).Utilizing empirical evidence from Li & Urai (2016), we implemented a rock salt rheology model that explains the interdependence between strain rate, stress, and temperature.Their data demonstrates that despite halite's appreciably higher effective viscosity compared to sulfur, it can accommodate dislocation creep.Notably, as temperature increases, the ductility of halite is enhanced, and effective viscosity decreases, aligning with the flow behavior proposed for Mercury's glacier-like formations.Specifically, in comparison to ice glaciers, which exhibit an effective viscosity of ∼10 12 Pa s as a result of dislocation creep leading to typical flow velocities ranging between 0.1 and 1 m day −1 , halite showcases distinct rheological characteristics.
Halite, subjected to ductile deformation at 200 °C, has an effective viscosity of ∼6 × 10 15 Pa s.This viscosity undergoes a significant reduction to ∼10 14 Pa s at 400 °C, similar to the viscosity of "polar type" (very cold) ice glaciers, thus signifying halite's marked thermal sensitivity.Consequently, even at 400 °C, halite sustains a viscosity roughly 100 times higher than that of temperate ice glaciers, suggesting potentially slower flow rates.Nevertheless, the exact runout distance as a function of time necessitates consideration of gravitational forces and the surface slope of the flow materials.
Given the marked rarity of Na 2 S in meteoritic and planetary materials, where its presence is confined mainly to caswellsilverite (NaCrS 2 ) and djerfisherite ) in enstatite meteorites and possibly on Mercury, our primary investigative focus pivots to NaCl.This compound is ubiquitously identified in various environments, including both aqueous and anhydrous settings, on Earth, Mars, and Ceres (Witze 2015) and has been identified in H-chondrite and ureilite meteorites (Berkley et al. 1978;Barber 1981;Zolensky et al. 1999).
Thus, we posit that NaCl, rather than Na 2 S, is a probable major component of the VRLs.The rheological behavior of halite aligns with the requirements for a primary compositional element in the generation of these glacier-like structures and possibly in the creation of hollows.However, this hypothesis does not eliminate the potential participation of other volatile substances.

Modeling Approach of Halite Flow
We base our model for runout from the initial emplacement of a ring mound on the Saint Venant equations (Miller 1984) for lubricated debris flow movement.The equations track vertically averaged mass conservation and momentum conservation along the horizontal axis.The governing equations are as follows: Mass conservation: with u(x, t = 0) = 0, h → 0 for x ?0; ∂h/∂x (x = 0, t) = 0, where h is the thickness of the mobile layer, u is horizontal velocity, x is distance, t is time, a = ∂z/∂x = local slope of the basement elevation z(x), and f is the angle of internal friction, which is approximately equal to the angle of repose for many situations.For many soils and rock types, the angle of repose ranges from about 25°-45°.Reduced gravity will lower the dynamic angle of repose.Following Kleinhans et al. (2011), under Martian gravity (which is almost the same as Mercury's), it has been estimated that large, very coarse-grained geomaterial has an angle of repose of about 31°.Our chosen value of 35°i s close to the Martian number for possible similar materials and is midrange for many other materials.The above mass and momentum conservation equations are in scaled form (that is, they are nondimensional), derived by first defining a characteristic length scale X o , a timescale T o , a depth scale H o , and a characteristic velocity scale U o .We then replaced unscaled variables (x, t, u, h) with the corresponding scaled variable (in italics) times the appropriate scale factor, e.g., x = x X o .Our scaled variables are: Values used for the characteristic scales are X o ∼ 10 4 m, H o ∼ 250 m, n o ranging from 10 3 -10 11 Pa s (for comparison, n H2O ∼ 10 −3 Pa s), and gravity g = 3.7 m s −2 .In the process of collecting terms in the scaling of the governing equations, we ended up with two nondimensional parameters, = e H X , and where e = ratio of depth to length scales, and b = ratio of viscous drag to gravitational pull.The initial condition is the profile of the uplifted crust and inward collapsed rim material.Additionally, the elevation z(x) of the base substrate is essential.The runout distance depends on the effective viscosity, gravity, and local slopes.Figures 3(D) and (E) include numerical solutions for two values of effective viscosity.The profile for lower viscosity (Figure 3(E)) tapers down to a slight elevation.The higher viscosity profile (Figure 3(D)) exhibits a blunted "nose" of the runout.In all cases, much of the profile is in the 100-250 m thickness range (similar to the profiles shown in Figure 2) after the volatile-rich material runs off the slope and eventually decreases to only tens of meters thick for the last kilometer of runout.Runout extends to about 10 km in these cases, also similar to profiles in Figures 2(B) and (C).Depending on the effective viscosity, the times corresponding to the late time profiles in Figures 3(D We posit that the Raditladi's formations, mirroring glacial characteristics, are a product of gravitational creep acting on exposed VRLs situated over a basaltic peak-ring topography (Figures 3(B) and (C)).A cornerstone of our findings lies in the rheological behavior of halite, which is consistent with the observed glacial deformation resulting from viscous flow.While we emphasize halite due to its abundance and widespread distribution in the solar system, as well as its ductile and volatile properties, it is pertinent to note, following Rodriguez et al. (2020) and Kargel (2013), that numerous other substances could play a part in Mercury's glacier-like flow formations and the emergence of hollows.Hence, halite serves as an illustrative but not exclusive model constituent, representing one of many potential contributors.
Halite as a glacial material (dynamics and emplacement timescales): Our computational model, considering a residual impact heat duration of approximately 1 Ma and a temperature range of 300 °C-500 °C (Figures 9(B), (C)), indicates that the halite VRLs would have experienced viscous flow due to the planet's gravitational influence (Figures 3(E), (D)).This process results in the formation of lobate profiles, with their characteristics being influenced by viscosity, as depicted in the halite-centric simulations (Figure 3(D)).
These features, exhibiting glacier-like behavior, consistently present runout distances predominantly ranging between 10 and 20 km (Figure 2(A)).This behavior aligns with our model's postulation that such propagation distances correlate with a thermal spike episode (Figure 9).The flow dynamics of halite also vary with the degree of differential stress applied.Under significant stress, halite demonstrates non-Newtonian flow characteristics.However, upon stress reduction, the flow behavior transitions to a Newtonian regime, marked by a linear relationship between stress and strain rate (Figure B2).
The flow dynamics of the material are initially accelerated due to the steepness of the peak ring slope, which prompts a significant reduction in effective viscosity.As this bulk material descends toward leveled plains, the driving force of gravity diminishes, leading to a slowdown in the creeping flow (Figures 3(C), (D); B2).Irrespective of the viscosity ranges considered in these simulations-from low to high-the thickness profiles are remarkably similar.They consistently show runouts extending to a distance of approximately 8-10 km, illustrating the consistent behavior of the material under varying conditions.
When considering the duration for the emplacement of glacier-like lobes, the internal temperature at the onset is a crucial factor.Our simulations suggest that if this initial temperature lies around 400 °C, within the predicted range of 300 °C-500 °C from iSALE2D simulations, then the process of deformation leading to the formation of these lobes would take roughly 150 Earth years (Figure B2).However, if the initial temperature is closer to 200 °C, reflecting the surface temperature at the conclusion of the iSALE2D simulations, the corresponding timescale extends to approximately 5000 Earth years (Figure B2).We note that these timescales are characteristic of large alpine glaciers found in both temperate and polar regions on Earth, reflecting similarities in the effective viscosities of halite on Mercury and ice on Earth.
Crucial observational findings constrain our age estimation.For example, we discerned the tectonic integrity of lobe margins overlaying faults on the basin floor (Figure 2).This superpositioning denotes a post-cooling and contraction phase placement of the lobe following the onset of the basin floor fracturing.In this context, we conjecture that Mercury's protracted nocturnal period of approximately 58 Earth days might have hastened the cooling process, instigating significant fracturing and paving the way for a potential millenniaspanning lobe progression.
The examined glacial surfaces lack tensional features, such as crevasses, which would be expected if active glacial dynamics had postdated the desiccation of the lithic layer.Furthermore, no discernible extensional tectonic patterning related to lobe spread is evident amid hollow distributions.
These findings imply that brittle extensional deformation within the lobes had largely ceased upon their surface desiccation, rendering faults improbable contributors to hollow formation.
Integrating these observations, we propose a model underpinned by a swift glacier-like advancement propelled by basal sliding and internal deformation processes.This active phase comes to a halt within a concise geological timeframe, possibly as brief as one to several Earth years.However, as our computational simulations suggest, it could have extended up to approximately 5000 Earth years.
From a broader perspective, our analysis conducted at the Raditladi site infers the presence of subterranean, volatileenriched stratigraphy within the upper crust extending beyond the demarcations of the planet's chaotic terrains documented in Rodriguez et al. (2020).This implication substantially expands the perceived extent of VRLs within Mercury's crust, indicating a planetary-scale distribution rather than constrained, localized occurrences.

Solar Heating Penetration Modeling
Solar irradiance on Mercury varies with latitude, daily cycles, and orbital phase, causing wide-ranging surface temperature fluctuations throughout a Mercurian day.Like Earth, where surface temperature variations affect only a few meters in depth, Mercury's subsurface temperatures likely remain stable daily and seasonally below a certain depth, despite a temperature gradient from geothermal energy flux.Thermal modeling can offer insights into these subsurface temperature profiles, even in inaccessible regions like Mercury.
To gauge the penetration depth of surface temperature variations into Mercurian regolith, we use a time-dependent 1D thermal conduction model, factoring in solar insolation changes over a Mercurian day (∼58 Earth days) as the top boundary condition.We target maximum penetration, focusing on equatorial conditions and neglecting latitudinal and longitudinal variations.We set the domain depth at 500 m, significantly deeper than the anticipated penetration depth.To ensure solution accuracy, we ran an energy-conservative finitedifference discretization of the domain at various grid resolutions.We incorporated geothermal heat flux at the model's base and surface radiation heat transfer.We solved for three thermal diffusivities, spanning an order of magnitude, to encompass the range of probable soil and rock properties.
The model equations and boundary conditions are included below.The initial condition is a steady-state solution with an average surface temperature (125 °C during one Mercurian day) and geothermal heat flux at the base.The temperature gradient can be described by ¶ where k T = K T /ρc p , k T is the thermal diffusivity, K T is the thermal conductivity, T is temperature, ρ is the bulk density, and c p is the specific heat.Values considered include K T = 3, 1.75 and 0.3 W m −1 K −1 , ρ = 3000 kg m −3 , and c p = 1000 J kg −1 K −1 .
At the surface, where the albedo α is 0.06, S(t) = S o sin(t/P) is the solar irradiance with a maximum irradiance, S o , of 9000 W m −2 , P is the rotational period of 58 Earth days, ε is the emissivityassumed to be 0.95, and σ is the Stefan-Boltzmann coefficient (5.6704 × 10 −8 W m −2 K −4 ).At depth (500 m here), where G is the geothermal heat flux, assumed to be 20 mW m −2 .Thermal properties internal to Mercury carry considerable uncertainty.Convective-conductive models of Mercury's thermal evolution (Breuer et al. 2007) suggest that the thermal conductivity of Mercury's mantle may be ∼4 W m −1 K −1 , but a lower value of 3 W m −1 K −1 may also be appropriate (Michel et al. 2013).Heat flow, combining conduction, radiogenic heating, and possibly convection yield a temperature gradient in the mantle of about 0.005 K m −1 (Figures 3 and 4; Breuer et al. 2007).If mantle thermal conductivity is as high as 4 W m −1 K −1 , then the heat flux to the surface would be 0.020 W m −2 , or 20 mW m −2 , the value used here in our simulation.Crustal thermal properties are expected to be different from mantle values.Given the uncertainties in crustal properties, we consider an order-ofmagnitude range (given above) in thermal conductivity.For a crust made of low porosity, competent rock, thermal conductivity will lie at the higher end of the range, while for a crust composed of high porosity, unconsolidated soil, thermal conductivity will be at the lower end.Further, the crustal density and specific heat are estimates.

The Potentiality of Prolonged Volatile Retention within Mercury's Shallow Subsurface
The Raditladi impact basin is thought to have formed ∼1 Ga (Prockter et al. 2009), and, according to our results, the development of its glacier-like features likely occurred within ∼1 Myr of the impact (Figures 9(B) and (C)).However, the hollows that have been identified in the region are believed to be relatively recent, perhaps currently active geological features, according to several studies (Blewett et al. 2011(Blewett et al. , 2013(Blewett et al. , 2016;;Helbert et al. 2013).
The preservation of the lobate shape of the glacier-like features, and their convex upward profiles in cross section with steep termini, as shown in Figure 2, suggests that most of the flow material must have been retained.Additionally, the presence of deep hollows indicates that a significant portion of these materials were volatile compounds, and their removal was localized to clusters.
We ran our simulations to ∼15,000 yr.By that time, the maximum depth of temperature variations stabilized (Figure 4(C)).As seen in profiles corresponding to the warmest and coldest parts of the day and night, the solar surface heating wave is characterized by a shallow penetration depth (∼3-8 m; depending on the geologic materials' thermal conductivity).Below that depth, there are only negligible diurnal fluctuations from the geothermal gradient (Figure 4(C)).Considering these shallow penetration depths, it is improbable that top-down temperature structural changes throughout the depth of the outflow layer would have significantly governed the dynamics of a fluidized bed.
This hypothesized thermal structure would account for the sustained bulk volatile retention evidenced by the hollows and moated craters.We posit that materials within the thermally affected depth produced a devolatilized upper layer, which effectively shielded and stabilized the underlying halite (and other volatiles) from solar heating effects (Figure 4(B)).
In the grand scheme of Mercury's geologic timeline, our findings imply that top-down devolatilization may have had a negligible role in depleting the volatile content of global surface and near-surface VRLs.This perspective suggests that the upper limits of the VRLs may have consistently remained shallow, a conjecture substantiated by observed discontinuities in impact crater rays and the presence of chaotic terrain hollows, as delineated in the study by Rodriguez et al. (2020).

Investigating Mercury's Polar Volatile-rich Stratigraphy:
Evidence and Implications The landscape of BC, showcased in Figure 1, is a fragmented, heavily cratered terrain marked by unique hummocky formations (Figure 8).These formations exhibit a scale-invariant pattern of crater diameters, with knobby rims that range approximately from 10-200 km (Figure 8).This pattern disputes the ejecta-origin theory for these terrains, as proposed by Ostrach et al. (2020).
Given the extensive regional burial, the ejecta process would have presumably wiped out smaller preexisting craters, and knobby rims should have been confined to the larger crater subset due to their burial by ejecta materials.However, as also observed within the chaotic terrain found antipodal to the Caloris basin (Rodriguez et al. 2020), the current state of BC, which includes widespread knobby craters of small diameter ranges, poses a challenge to this theory.
We find more merit in an alternative hypothesis that volatilerich constituents within the upper crust underwent sublimation, causing fragmentation and disaggregation of the regional landscape.This interpretation is consistent with previously posited mechanisms for the chaotic terrain found antipodal to the Caloris basin (Rodriguez et al. 2020).
The base of disaggregation within Borealis Chaos showcases uniformly leveled plains, implying a region-wide uniformity in maximum depth of relief loss.We posit that these plains demarcate the stratigraphic contact between the gravimetrically detected possible lava-covered cratered paleosurface (Figure 8), previously documented by Qingyun et al. (2018), and the regional base of a VLR.In other words, our interpretation situates this paleosurface as the upper boundary of a volatiledevoid basement (Figure 10(A)).In Qingyun et al.ʼs (2018) study, the gravity signatures' source ambiguity was narrowed by limiting the bandwidth to 12-75, pointing to a possible crustal origin.Based on this constraint and the ring-shaped anomalies, they suggested these were ancient impact structures.Furthermore, the spatial density analysis of the gravity anomalies and ghost craters allude to an epoch predating the LHB, indicating the presence of a buried cratered paleosurface, potentially underlain by lavas.

Deciphering the Genesis and Evolution of Mercury's
Volatile-rich Layers: Evaluating Potential Atmospheric Sources and Mechanisms We hypothesize that the proposed cratered paleosurfaces (Qingyun et al. 2018) concealed beneath Mercury's VRLs (Figure 10(A)) demarcate a transition phase, wherein volatile emplacement followed the solidification and crater retention capacity of the planetary surface.However, this hypothesis is inconsistent with the VRLs originating from a buoyant sulfide primary crust formed through the solidification of a sulfide-rich magma ocean (Parman et al. 2016).Hence, we propose an alternate genesis, characterizing it as a mega-condensation feature predating the LHB.We further suggest an origin by atmospheric collapse, a process recognized for generating global volatile layers on planetary surfaces (Soto et al. 2015;Sullivan 2015;Arimatsu et al. 2020).
Our results imply the complete disaggregation of surfaces that are considered to have experienced cratering during and after the LHB (Ostrach et al. 2020).Consequently, we infer that the VRLs are vestiges of an archaic stratigraphic regime antedating the LHB epoch.Atmospheric collapse (virtually complete condensation) of volatiles is a possible mechanism.The volatiles could have included ice-forming materials such as H 2 O, slightly volatile materials such as NaCl, and many others.The stratigraphy, embedded between the formation of the gravity-deduced lava-shielded crater surface and the LHB, indicated that the atmospheric collapse must have occurred within a narrow geologic timeframe.Notably, had there been an extended, slow buildup of VRL materials, this likely would have led to alternating stratified zones of volatile and nonvolatile (presumably volcanic) materials.Craters created within this layered formation would have included subsets with VRL-composed rims and others with basaltic rims; hence, the disintegration of heavily cratered terrains would have affected only the VRL subsets, a circumstance we did not observe.
We propose the rapid development of a Mercurian atmosphere (Jäggi et al. 2021), potentially driven by large Kuiper Belt object impacts (Bekaert et al. 2020), volcanism (Kasting 1993), or impact-induced excavation into a gas-charged mantle (Marchi et al. 2016).The Kuiper Belt object hypothesis is consistent with a proposed phase of pre-LHB asteroidal bombardment (Hartmann 1984(Hartmann , 2019) ) and the finding of a possible cratered paleosurface buried beneath the heavily cratered terrains (Qingyun et al. 2018).
In addition, this early bombardment likely engendered an early Mercury megaregolith analogous to the lunar counterpart (Liu et al. 2022).Consequently, we propose a novel potential atmospheric source-the sudden degasification of this primordial megaregolith.Given the expected high porosity of this layer, it is scientifically plausible to consider it a significant reservoir for volatiles.Drawing geological parallels, the destabilization of a similar megaregolith on Mars may have contributed to the genesis of a northern ocean approximately 3.4 Ga (Clifford 1993;Clifford & Parker 2001).Furthermore, Mercury and Mars display similar, enriched K/Th and Cl/Th ratios compared to Earth and Venus, hinting at parallel trends in their volatile content evolution (Nittler & Weider 2019).
We propose that the swift degassing, hypothesized to have contributed to Mercury's primordial atmosphere (Jäggi et al. 2021), was likely synchronous with accelerated condensation during prolonged night periods.The frigid nocturnal conditions probably incited atmospheric instability, catalyzing the condensation and emplacement of volatile compounds; hence, providing a mechanism to lead to the development of the VRLs.Following volcanic events (Marchi et al. 2013), the VRLs may have been intermittently buried in areas, accounting for areas where chaotic terrains are less prevalent (Figures 10(B)-(D)).Within this theoretical construct, the glacier-like structures observed on the Raditladi peak ring (Figure 2) could represent zones where VRL constituents have been unearthed (Figures 10(B)-(D)).

Possible Origins of Halite-Dominated Volatile-rich Layers
The genesis of halogen-rich minerals in planetary crusts, including halite, is commonly attributed to the aqueous alteration of Na and Cl-containing minerals, such as select feldspathoids, micas, and amphiboles (Volfinger et al. 1985;Jenkins 2019).Alternatively, these minerals can also form through primordial and subsequent HCl and NaCl magmatic exhalations (Harlov & Aranovich 2018;Ueda & Shibuya 2021).The role of water in volcanic exhalations can be significant, as observed, for example, in erupted materials from Mount Saint Helens and other volcanoes, which include halite-and other salt-encrusted ash particles (Varekamp et al. 1986;Staunton-Sykes et al. 2021).
Earth's mantle halogens are predominantly released through volcanic emissions, which include gases, liquids, and chalcogens, with metals present in trace amounts (Millard et al. 2006;Guo & Korenaga 2021).Consequently, vapor redistribution of exhalative salts from volcanic activity, a phenomenon possibly particularly extensive and intensive on Mercury due to its comparatively hotter surface conditions, could have contributed to VRL formation.
For example, Garrison et al. (2000) investigated the chlorine abundance in meteorites, specifically focusing on 14 enstatite chondrites.These meteorites, considered as likely analogs to Mercury's precursor material according to Rodriguez et al. (2020), revealed a chlorine concentration varying between 24 and 653 ppm, with a mean value of 257 ppm.Notably, this average markedly surpasses the chlorine concentration found in the bulk silicate of Earth, which is reported to be 36 ppm (Kargel & Lewis 1993).
Upon investigating the mean concentration of Cl found within enstatite chondrites and assuming stoichiometric conversion of Cl and Na into NaCl, we can extrapolate potential NaCl abundance.In this stoichiometric transformation, we assume the direct conversion of available Cl and Na in the source material into NaCl, predicated on their relative atomic masses and the theoretical 1:1 molar ratio of Na:Cl in NaCl.The actual abundance of NaCl, however, could deviate from this calculated maximum due to the heterogeneous distribution of these elements and varying reactivity toward other potential cationic species in the system.
The projection of these calculations onto Mercury's global scale-factoring in the planet's mass, surface area, and overall Cl content, we find that a complete outgassing and uniform condensation of NaCl could hypothetically yield a global layer of approximately 857 m thick.Although idealized, this estimation presents an order-of-magnitude approximation supported by the potential outgassing of other volatile species.In addition, Piani et al. (2020) provided evidence of significant water content in enstatite chondrites, suggesting a potential contribution to the formation of large water bodies.This finding suggests that an estimated 857 m thick, global VRL could represent a mega-precipitation feature from transiently ponded brines, potentially forming extensive lakes with an extended stability under the faint young sun.These brines could have been released through volcanic-hydrothermal activity, Kuiper Belt object impacts, or even concurrently with the rapid degasification of the primordial megaregolith.
Details on Theoretical Estimation of NaCl Layer Thickness on Mercury.Mercury's surface, known as regolith, is composed of a mixture of fine dust and small rocks and contains a mean chlorine concentration of 257 ppm.Considering that all this chlorine could react with sodium to form NaCl, we could envisage a potential NaCl layer on Mercury's regolith.The estimated thickness of this theoretical NaCl layer would be approximately 857.1 m. Procedure: The hypothetical thickness of this NaCl layer is computed via the following equation: The results of this theoretical calculation, which presupposes all chlorine on Mercury as available to form NaCl, suggest a potential NaCl mass of approximately 2.57 × 10 14 grams, yielding a consequent global mean layer thickness of about 857.1 m.
Limitations: Our model posits chlorine's existence in a minimally complex combined state as NaCl; however, before forming NaCl, Cl would have been in other mineral forms and might not have been entirely available to form NaCl. Furthermore, the interplay of released Cl − or Cl 2 with diverse elemental constituents could significantly modulate the actual thickness, composition, viscosity, and volatility of the hypothesized NaCl-rich stratum on Mercury's regolith.

Geo-energetic Forces in Mercury's Polar Disaggregation
The observation of extensive north-polar plains has been interpreted as a series of flood basalt-like volcanic occurrences after the LHB (Head et al. 2011).An increased geothermal flux linked with this phase of volcanic activity could have destabilized north-polar VRL deposits, triggering the formation of BC.However, this hypothesis is challenged by the absence of supporting evidence for north-polar volcanism, such as discernible vents (Ostrach et al. 2020).This incongruity compels us to consider alternative geo-energetic factors that could have induced VRL disaggregation, such as progressive solar heating (Feulner 2012), which has been previously considered as a mechanism for chaotic terrain formation, potentially leading to continuous disaggregation over several billion years (Rodriguez et al. 2020).
Recent research shows that magnetic cusps in Mercury's north polar region could be sputtering north-polar volatile compounds (Raines et al. 2022), suggesting the potential of ongoing VRL erosion coupled with high-energy particles.Furthermore, during its early geological phase, Mercury is believed to have possessed a strong geomagnetic field, similar in strength to Earth's (Johnson et al. 2015), and we suggest that its configuration could have led to high rates of regional VRL sputtering due to focused high-energy particle influx.
For example, Mercury's magnetic field north-south asymmetry has its equator northward of the geographic equator, a configuration that could have amplified solar wind interactions, escalating the rate of magnetic reconnection events over the northern hemisphere.These events, probably coupled with super solar flares from the primordial Sun (Killen et al. 2012), could have precipitated charged particles along north-polar converging magnetic field lines, focusing sputtering and erosion.This hypothetical scenario explains why chaotic terrains are absent from the south pole.In this scenario, nonpolar chaotic terrains could have developed in connection to geomagnetic storms instigated by solar phenomena such as coronal mass ejections and solar flares, which might have disrupted Mercury's magnetic field lines, potentially generating regions of high-energy particle bombardment due to inadequate geomagnetic shielding.
The proposed association between solar particle events and the formation of Mercurian chaotic terrains explains the heterogeneous distribution of these terrains, which deviates from purely regional geological contexts.Persistent magnetic disruptions, whether due to local anomalies, solar wind funneling, or magnetospheric 'cusps,' are likely to have played a significant role in chaotic terrain formation.It is important, however, to consider that these events could have initially increased surface shortwave roughness, potentially accelerating erosion due to uneven solar heating compared to more topographically homogenous surfaces.In addition, Mercury's chaotic terrains, characterized by hummocky landscapes, suggest a volume reduction of up to a few kilometers.If the regional VRL is dominated by halite, incoming solar protons could instigate radiolysis, leading to its dissociation into ions or radicals such as Na and Cl.This scenario provides an effective erosion mechanism congruent with the high abundance of these elements observed in the exosphere.

Implications toward Cosmoglaciology and Habitability
Our discovery of potential glaciation on Mercury, along with the recent identification of nitrogen glaciers on Pluto (Stern et al. 2015), expands our understanding of glacial phenomena beyond the familiar water-ice glaciers of Earth and Mars, reaching into both the outermost and innermost regions of our solar system.Notably, the likely occurrence of salt glaciers on Mercury, Earth, and Mars suggests that glacier-forming salts may be highly adaptable to various environmental conditions.
Research conducted in Earth's Atacama Desert lends credence to the proposition that life can persist within hydrated salt nodules, forming ecological enclaves that mimic the planetary environmental Goldilocks zones capable of sustaining liquid water (Davila et al. 2010(Davila et al. , 2015;;Davila & Schulze-Makuch 2016).Similarly, on Mars, substantial hydrological discharges have been causally connected to the thermal dehydration of salt deposits (Kargel et al. 2007).These findings have significant implications for the search for life in extraterrestrial environments, especially those featuring similar geochemical landscapes.
Paralleling these Earth and Mars examples, Mercury's VRLs could present a comparable scenario where zones of apt temperature and pressure host stable niches of saline water, potentially habitable environments.Intriguingly, impact-heated halite on Mercury may have safeguarded these saline lenses for an extended duration, approximately ∼1 Myr.A thorough interpretation of the temperature gradient across a ∼2 km thick VRL is pivotal to assessing the possibility of the sustained existence of liquid water on Mercury, thereby enriching our perception of these subsurface Goldilocks zones.
Our study underscores the significance of subsurface cryospheres in terms of their age and extent.For instance, we find that Mercury's VRLs, extending ∼2 km in thickness near the equator (Rodriguez et al. 2020), likely formed before the LHB (>3.8 Ga).Similarly, Mars' cryosphere, with an equatorial thickness ranging from ∼2.3-4.7 km (Clifford & Parker 2001), is probably associated with global volatile cycling ∼3.7 Ga or earlier (Clifford & Parker 2001;Rodriguez et al. 2015).The inferred extremely ancient ages of these cryospheres, aligning with the timeline of life's emergence on Earth, coupled with their extensive spatial scales, suggests the potential for analogous evolutionary processes and biological diversification, albeit within a subsurface context.In addition, the large vertical and horizontal domains suggest that hydrothermal circulation could have enhanced the sustainability of habitable conditions by continuously supplying fresh brine and minerals (Travis et al. 2013).
Finally, the growing number of exoplanets discovered in close orbits around their stars suggests that cryospheres similar to Mercury's could be a common feature in our galaxy.Therefore, we suggest that glaciation could be a prevalent phenomenon on exoplanets with unstable cryospheres, offering new perspectives on the complexity and diversity of planetary environments and establishing the prospect of common subsurface Goldilocks zones.

Summary and Implications
Detecting widespread elemental volatile surface compositions, ubiquitous sublimation hollows, and extensive chaotic terrains has significantly reshaped our perception of Mercury's geological past.These observations collectively point to the presence of volatile-rich strata spanning several kilometers in depth and likely formed before the LHB (∼3.8 Ga).This notion challenges the conventional view of a volatile-depleted Mercurian crust.
The morphologies within Mercury's Raditladi basin bear a striking morphologic resemblance to glaciers on Earth and Mars, suggesting their origin from an impact-exposed, VLR, likely containing halite.Our numerical simulations show that the unique rheological properties of halite, including the high thermal sensitivity of its viscosity, reinforce this hypothesis.These glacier-like features occur beyond the chaotic terrain boundaries, indicating a potentially global yet concealed, volatile-rich upper stratigraphy.We posit that the exposure of these volatile-rich materials, instigated by impact events, could have been instrumental in the formation and evolution of hollow features, signifying a complex geodynamic history of volatile migration and redistribution, essentially interconnecting some of the oldest and youngest stratigraphic materials on the planet.
Our analysis of Borealis Chaos indicates a regional extension of a VRL over the north pole, providing unique stratigraphic insights into its genesis.We observe a regionally uniform collapse base, interpreted as a likely transition to a volatile-free, presumably volcanic basement.This basement, evident in previous gravimetric data set analyses, is likely an impacted, lava-covered paleosurface.The VRL's superposition on this cratered terrain necessitates a deposition-based megastructure hypothesis, potentially resulting from the collapse of a transient, primordial atmosphere.This early atmosphere could have originated from exogenic sources like comets or endogenic processes such as volcanic degassing or sudden degasification of a primordial megaregolith.Given the VRLs' possible halite-dominated composition, we also consider an origin due to a mega-precipitation feature from extensively ponded brines (e.g., transient salty lakes, seas).
Three potential mechanisms are proposed for the disintegration of VRLs in Mercury's northern polar region: post-LHB magmatic intrusion, combined impact of solar flares and extreme solar events with a stronger dipole magnetic field, and gradual intensification of solar heating.In addition, extensive volcanic activity could have effectively concealed these VRLs throughout other regions of the planet, providing a possible explanation for the relatively sparse distribution of chaotic terrains on Mercury and their absence in the southern polar region.Hence, our assessment is that Mercury's VRLs were initially emplaced as a global stratigraphic deposit, which might still retain much lateral continuity.
We assigned values for surface height and albedo to each maplet vertex.We used stereo positioning to determine the height of the maplet's center vertex.We identified a corresponding pixel in each image corresponding to the exact physical location on the surface.Then, we extended rays from each pixel and selected the position that minimized the distance error to each ray, with higher-resolution images weighted more heavily in the process.The cluster where the rays converged identified the center vertex position of the maplet.We used photoclinometry to define the height and albedo of the remaining maplet vertices from the images.The brightness of each image pixel was used to determine the albedo and slope of the surface.A bright pixel could be due to high albedo or surface tilting toward the spacecraft camera.A dim pixel could be due to low albedo or surface tilting away from the camera.We found that high-incidence images emphasized the albedo of the surface, while low-incidence images revealed the surface's topography.Using photoclinometry, we solved for the albedo and slope for each maplet vertex by using images with different incidence angles.We then integrated the height of the center maplet vertex and the slope of all other maplet vertices to determine the height for all remaining maplet vertices.
We used hundreds of maplets to cover the DTMs created for this study, with each maplet having a portion of topography that overlaps with adjacent maplets.We used the heights of overlapping maplets to ensure that the photoclinometric uncertainty remained small at the edges of each maplet.These overlaps are crucial, as any errors in the slope will increase the uncertainty as the height solution progresses outward from the stereo-defined center vertex.
After generating the maplets, the following necessary process in SPC is to update the image's spacecraft position.Once the position and height of many maplets had stabilized, we performed a reverse of the stereo solution.Instead of using multiple images to define the center position of a maplet, we used many maplets to define the image's spacecraft position.We then used the updated spacecraft position to refine further the stereo position of each maplet, which in turn was used to refine the spacecraft position.Thus, SPC performs an iterative solution between the topography and the spacecraft position.During all of these iterative processes, we continually improved the alignment of images within each maplet.When the topography and spacecraft position no longer improved with successive iterations, we combined the maplets into a single DTM for distribution.Further SPC details can be found in Gaskell et al. (2008Gaskell et al. ( , 2023)), Palmer et al. (2022).SPC error analysis can be found in Weirich et al. (2022).Details of the three DTMs are shown in Table A1, with visual representations in Figures A7-A9    The profiles illustrate the interrelation between flow thickness and runout length at discrete temporal junctures.The corresponding times (t i ) depend on the initial internal temperature of the flowing layer.Given an initial internal temperature approximating 400 °C, falling within the 300 °C-500 °C range inferred from the iSALE2D simulations, the timescale (t 4 ) is estimated to be of the order of 150 yr.Alternatively, if the initial internal temperature approximates 200 °C, mirroring the surface temperature at the conclusion of the iSALE2D simulations, the timescale (t 4 ) is projected to span roughly 5000 yr.Upon the flowing material entirely vacating the peak ring slope, any ensuing creep flow is anticipated to markedly decelerate.

Figure 2 .
Figure 2. (A) Morphological map of the Raditladi impact basin (27°.28N, 240°.93W), highlighting the distribution of hollows (yellow), lobate debris aprons bordering the peak ring ridges (red), fractures within the plains enclosed by the peak ring (white lines), and small, moated craters (blue dots).(B) and (C) Detailed views of a lobate feature, with yellow arrows marking its frontal margin and red arrows pointing out clusters of hollows.Panels (D) and (E) illustrate elevation profiles X-X′ and Y-Y′, respectively, indicating the lobe's projected maximal thicknesses of ∼200 and ∼250 m.These estimates are based on the assumption that the nearly flat terrain at the lobe's margin (red-circled white dot) persists under the lobate deposit.The pre-hollow topography is approximated by linking the surrounding elevations.This reconstruction suggests the hollows generated a negative relief of about 50-75 m (areas marked by red lines in profile Y-Y′).Data sources: Panel (A) consists of a darkened MESSENGER MDIS global base map (∼166 m px −1 ; courtesy of NASA/JPL/GSFC).Panel (B) integrates a geo-registered MESSENGER Regional Targeted Mosaic (RTM N01_003828_2238046) and a custom DTM (DTM RAA201_EQv1), whereas panel (C) combines geo-registered RTMs (N01_003828_2238046 & N01_000167_0593484) with a custom DTM (DTM RAB202_EQv1).Detailed information on RTM processing and DTM creation can be found in the supplementary materials.

Figure 3 .
Figure 3. (A) Peak ring sections with flat tops, marked by red arrows.Elevation profiles L-L′ and Z-Z′ chart the topographic changes over steep-sided and flat-topped peak ring sections, respectively.We utilized their topography to generate the profiles for the peak ring's VRL-shed sequence illustrated in panels (B) and (C).(B) The schematic diagram illustrates the reconstructed stratigraphy of the peak ring, including a VRL zone atop the impact-uplifted basaltic crust.Drawing from our numerical data in Figures 9(B) and (C), we suggest that the impact heat induced a thermal anomaly beneath the peak ring, leading to a partially melted region (medium-tone orange).VRL materials with minimal melt are depicted in light-tone orange.(C) Diagram correlating flat-topped mesas with the hypothesized source of morphologic glaciers.Our theory suggests instability in the peak ring ridges' upper sections caused ductile flow, leading to the creation of expansive lobate aprons that share morphological traits with glaciers (dark-tone orange).(D), (E) Flow simulations account for factors such as stratigraphic composition, thermal structure, topography, and volatile-rich composition.The simulations display lobate profiles in gravity-induced flow, with runout distance over time being dependent on variables like viscosity, gravity, and substrate slope.The simulation reveals an elevation range of 100-250 m, tapering down to a few tens of meters for the S 8 case (E) and a steep front for the halite case (D).The labels T 1 − T 5 denote times of flow evolution, with T 0 representing the initial thickness of the VRL.Late-stage profiles range from hours to days for the S 8 flow, while they stretch to millennia for the halite flow.Our simulated profiles closely mirror the observed ones depicted in Figure 2(B).Data sources: (A) These reconstructions incorporate geo-registered RTMs (RTMs N01_000167_0593484 & N01_003828_2238046) and a DTM (DTM RAD202_EQv1).Comprehensive information about the RTM processing and DTM generation can be referenced in the supplementary materials.

Figure 4 .
Figure 4. (A) Moated craters (white arrows) on a morphologic glacier-like lobe, as captured in the geo-registeredRTM (image N01_003818_2236022). Detailed procedures for RTM processing can be found in the supplementary materials, with context and geographic location provided in Figure 2(A).(B) Cross-sectional reconstruction, drawn along line M to M′, postulates the formation of moats due to impact-induced exposure of near-surface, volatile-rich materials along the edges of some small (subkilometer to a few kilometers) craters.(C) The modeling of surface solar heat waves reveals a relatively shallow penetration depth (∼3-8 m), consistent with the long-term retention of volatiles (see Sections 4.2.2.4 and 4.2.2.5).The presented profiles represent the thermal extremes, with the warmest portion of daytime (upper curves) and the coldest point of nighttime (lower curves).This finding, aligning with the hypothesis that hollows are recent geological features, extends the same process to explain the formation of moated craters.The model incorporates a range of parameters, including a peak solar irradiance of 9000 W m −2 , emissivity set at 0.95, an albedo of 0.06, and a geothermal flux of 20 mW m −2 .Thermal conductivities are specified as K t = 3 W m −1 K −1 for intact rock, K t = 1.75 W m −1 K −1 for a composite substance, and K t = 0.3 W m −1 K −1 for highly porous soil.

Figure 6 .
Figure 6.(A) Eminescu crater's central peak surrounded by high albedo lobate structures.(B) The enlarged view reveals the bright regions as clustered depressions with marginal patchy textures.We suggest these depressions originated from the loss of volatiles in the VRL lobes extending from the peak ring.Data sourced from MESSENGER's QuickMap tool Regional Targeted Mosaic (RTM).

Figure 5 .
Figure 5. (A) Overview of the Raditladi impact basin, indicating the region magnified in panel (B), which zeros in on a landslide extending from the southeastern rim (highlighted by a white arrow).Notably, the landslide lacks associated hollows.For comparison, panel (C) introduces the Holbeck landslide from North Yorkshire, England-a coastal landslide-marked by a white arrow.Data source: panel (A) consists of a MESSENGER Mercury global DEM (∼665 m px −1 ) draped over a MESSENGER MDIS global base map (∼166 m px −1 ; courtesy of NASA/JPL/GSFC).Panel (B) consists of an MDIS base map (∼166 m pixel −1 ; courtesy of NASA/JPL/GSFC).

Figure 7 .
Figure 7.View showing the locations of north polar high gravity anomalies interpreted to be locations of craters and basins that were covered by thick stacks of northern-plain-forming lavas before the LHB (Qingyun et al. 2018).The color-coded, band-limited, degree-strength Bouguer anomaly map of the northern high latitudes of Mercury (regions above 45°N) is superimposed on the MESSENGER MDIS global base map.The MDIS base map furnishes a higher spatial resolution, approximating 166 m per pixel.

Figure 8 .
Figure 8. (A) View portraying the BC landscape, characterized by its segmentation of heavily cratered terrains.Notice the similarly degraded craters of varying diameters (10-160 km) exhibiting rim discontinuities (arrows 1-3).The dashed curve traces the location of the profile Z-Z′.(B) A close-up of BC's floor reveals a crater-saturated surface, including notably degraded craters morphing into irregular depressions (feature 1).Larger irregular depressions (feature 2) and even more extensive plains with enclosed mesas consisting of heavily cratered plateau residues (feature 3) can also be seen.A regionally consistent topographic base around −4000 m indicates the interface between a regional VRL and the gravity-detected cratered, lava-covered paleosurface (refer to Figure7).Data source: both panels consist of a MESSENGER Mercury global DEM (∼665 m px −1 ) draped over a MESSENGER MDIS global base map (∼166 m px −1 ; courtesy of NASA/ JPL/GSFC).
) and (E) range from several hours to days for an S 8 -dominated flow or many years for a halitedominated flow.Viscous creeping flows on Earth provide examples of relevant flow rates(Takagi 2010).Lava flow and rock slide rates on Earth can range from several meters per second to less than 1 m day −1 on small slopes(Takagi 2010).The effective viscosity of various mafic lava flows on Earth is in the range of 10 3 -10 6 Pa s(Chevrel et al. 2019; water has a viscosity of ∼10 −3 Pa s at 20 °C).Ice glaciers deform by slow creep of dislocations and exhibit an effective viscosity of roughly 10 12 Pa s with typical flow rates of 0.1-1 m day −1 .The effective viscosity of ductile halite flow at ∼200 °C (Li & Urai 2016) is about 6 × 10 15 Pa s.At 400 °C, halite viscosity drops to about 10 14 Pa s. 4.2.2.3.Modeling Results of Glacier-like Features: Rationale for Halite Glaciation and the Possibility of Global Volatile-rich Upper Crustal Stratigraphy

Figure 10 .
Figure 10.Illustrative diagrams depicting key geological scenarios: (A) Gravity data (Qingyun et al. 2018) reveals the presence of an ancient cratered terrain, which we propose to have been overlain by a global VRL.This panel shows a cross-sectional view of this stratigraphy along a longitudinal reconstruction spanning from the equator to the north pole.(B) Extensive areas of the VRL are concealed beneath basalts, while the VRL remains exposed in the north polar region.(C) Impact events on the basaltic surface result in crater formation and the exhumation of buried VRL materials.Craters are also observed in the north polar area.(D) The exhumed VRL material flows and gives rise to glacier-like lobate features.The north polar region undergoes collapse, leading to the degradation or elimination of most craters.
Thickness = (Mass of NaCl)/(Area of Mercury × Density of NaCl) Where 1. Mass of NaCl: This term represents the potential mass of sodium chloride that could form, calculated by multiplying Mercury's mass by the chlorine concentration and then by the ratio of the molar mass of NaCl to that of chlorine.2. Area of Mercury: The total surface area of Mercury amounts to 7.48 × 10 11 m 2 .3. Density of NaCl: Sodium chloride is 2.17 g cm −3 , or equivalently, 2170 kg m −3 .The Mass of NaCl is elaborated as: Mass of NaCl = (Mass of Mercury) * (chlorine concentration) * (molar mass of NaCl/molar mass of chlorine) On substituting known values, we derive: Mass of NaCl = (3.301× 10 23 kg) * (257 ppm) * (58.443 g mol −1 /35.453 g mol −1 ) = 2.57 × 10 14 g Results and Implications:

Figure A1 .
Figure A1.The image depicted above, taken by NASA's orbital MESSENGER spacecraft, underscores the faint blue shade (here in exaggerated color contrast) of unique geological formations termed "hollows" on Mercury.Mercury's landscape is punctuated by impact craters consisting of bowl-shaped depressions, each typically accentuated by raised rims and accompanying ejecta blankets.In stark contrast, hollows emerge as irregular, pit-like structures, their steep walls descending into flat or gently sloping floors.An unmistakable feature of these hollows is the surrounding bright halos, a characteristic that sets them distinctly apart from craters, as they notably lack the common traits of raised rims and ejecta blankets.The showcased image offers a detailed view of a 40 km segment of the Raditladi impact basin's base, featuring its striking central peak mountains.Image Credit: NASA/JHU APL/CIW.

Figure A2 .
Figure A2.Comparative Perspectives of Different Glacier Types-Raditlati's lobe presumed to be a halite glacier (A), a debris-covered Martian ice glacier (B) at coordinates 102°58′E, 40°36′S, and a salt glacier in Iran's Zagros Mountains (C) at coordinates 54°54′E, 27°59′N.The white and red arrows distinctly point toward the originating mountains and their respective lobate glacier fronts.

Figure A3 .
Figure A3.Sequential representation of iSALE2D hydrocode model results.(a) Snapshot at 30 s post-impact, illustrating the maximum depth achieved.(b) Frame at 170 s, depicting the ejecta curtain encapsulating the overturned stratigraphy.(c) Image at 300 s showcasing the formation of a transient central peak and ejecta blanket.(d) Final visualization at 700 s (see also Figure 9(A)).

Figure A4 .
Figure A4.iSALE2D hydrocode model incorporating a dunite crust-A potential superior match to basalt, assuming that Mercury's crust shares closer thermal and physical characteristics with komatiites.Observe that the peak ring area, approximately 60 km from the center, continues to harbor VRL materials.(a) Snapshot at 30 s post-impact, illustrating the maximum depth achieved.(b) Frame at 170 s, depicting the ejecta curtain encapsulating the overturned stratigraphy.(c) Image at 300 s showcasing the formation of a transient central peak and ejecta blanket.(d) Final visualization at 700 s (see also Figure 9(A)).

Figure A5 .
Figure A5.iSALE2D hydrocode model incorporating an initial regolith overburden-The model features an overlying stratum of 1.2 km of basaltic regolith placed atop the VRL.Take note of the peak ring, approximately 60 km radially distant from the center, which continues to preserve VRL constituents.(a) Snapshot at 30 s post-impact, illustrating the maximum depth achieved.(b) Frame at 170 s, depicting the ejecta curtain encapsulating the overturned stratigraphy.(c) Image at 300 s showcasing the formation of a transient central peak and ejecta blanket.(d) Final visualization at 700 s (see also Figure 9(A)).

Figure A6 .
Figure A6.Implementation of iSALE2D hydrocode model with enhanced initial VRL-This model features a considerably thickened VRLr, extending to a depth of 11 km.Note the preservation of VRL constituents within the peak ring, located approximately 60 km from the central point.(a) Snapshot at 30 s post-impact, illustrating the maximum depth achieved.(b) Frame at 170 s, depicting the ejecta curtain encapsulating the overturned stratigraphy.(c) Image at 300 s showcasing the formation of a transient central peak and ejecta blanket.(d) Final visualization at 700 s (see also Figure 9(A)).

Figure A7 .
Figure A7.(A) MESSENGER MDIS Image N1015454462M showcasing the DTM location centered at coordinates 28°.340 Lat, 120°.283 E Lon. (B) A DTM rendering using shaded relief or hillshade, providing elevation information.(C) DTM rendering under the identical incidence and emission angle as the spacecraft image, delivering both elevation and albedo data.All panels in the figure span an area of 17 km across. .

Figure A8 .
Figure A8.(A) MESSENGER MDIS Image N1046021896M illustrating the DTM location centered at coordinates 27°.620 Lat, 120°.442 E Lon. (B) A hillshade rendition of the DTM, solely presenting elevation data.(C) A DTM rendition that matches the incidence and emission angle of the spacecraft image, thereby offering both elevation and albedo details.Each panel within this figure extends over a span of 17 km.

Figure A9 .
Figure A9.(A) MESSENGER MDIS Image N1045992946M showcasing the DTM location centered at coordinates 27°.944 Lat, 120°.277 E Lon. (B) A hillshade interpretation of the DTM, exclusively providing elevation details. (C) DTM visualization conducted under the same incidence and emission angle as the spacecraft image, revealing both height and albedo aspects.Each panel within the figure spans an area of 17 km across.

Figure B2 .
FigureB2.The image delineates the varied behavior of halite under differing viscosity conditions during simulations.When the viscosity is elevated, the resultant flow exhibits a sharp leading edge, as depicted in the figure.The profiles illustrate the interrelation between flow thickness and runout length at discrete temporal junctures.The corresponding times (t i ) depend on the initial internal temperature of the flowing layer.Given an initial internal temperature approximating 400 °C, falling within the 300 °C-500 °C range inferred from the iSALE2D simulations, the timescale (t 4 ) is estimated to be of the order of 150 yr.Alternatively, if the initial internal temperature approximates 200 °C, mirroring the surface temperature at the conclusion of the iSALE2D simulations, the timescale (t 4 ) is projected to span roughly 5000 yr.Upon the flowing material entirely vacating the peak ring slope, any ensuing creep flow is anticipated to markedly decelerate.