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

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):

$$\begin{eqnarray*}&&{D}_{r}=0.91\displaystyle \frac{{D}_{\mathrm{tr}}^{1.125}}{{D}_{Q}^{0.09}}.\end{eqnarray*}$$

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⁻¹ 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).

Zoom In Zoom Out Reset image size

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, r–z 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, 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).

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)).

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₂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₈) and NaCl as the likely predominant constituents of the glacier-like features. Here, S₈ 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₂. However, these ranges intersect with the melting point spectrum of sulfur (S₈), 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₈ (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¹² Pa s as a result of dislocation creep leading to typical flow velocities ranging between 0.1 and 1 m day⁻¹, halite showcases distinct rheological characteristics.

Halite, subjected to ductile deformation at 200 °C, has an effective viscosity of ∼6 × 10¹⁵ Pa s. This viscosity undergoes a significant reduction to ∼10¹⁴ 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₂S in meteoritic and planetary materials, where its presence is confined mainly to caswellsilverite (NaCrS₂) and djerfisherite ${({{\rm{K}}}_{6}\mathrm{Na}({\mathrm{Fe}}^{2+})}_{24}{{\rm{S}}}_{26}\mathrm{Cl})$ 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₂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.

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:

$$\begin{eqnarray*}&&\partial h/\partial t+\partial ({uh})/\partial x=0,\ h(x,t=0)={h}_{o}(x)\end{eqnarray*}$$

Momentum conservation:

$$\begin{eqnarray*}\begin{array}{rcl}\partial u/\partial t+u\partial u/\partial x & = & \sin (a)-e\partial h/\partial x+b{\partial }^{2}h/(\partial {x}^{2})\\ & & -\,\cos (a)\tan (f)\end{array}\end{eqnarray*}$$

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° is 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:

$$\begin{eqnarray*}&&x=x/{X}_{o},\end{eqnarray*}$$

$$\begin{eqnarray*}&&h=h/{H}_{o},\end{eqnarray*}$$

$$\begin{eqnarray*}&&u=u/{U}_{o},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{U}_{o}={X}_{o}/{T}_{o},\end{eqnarray*}$$

$$\begin{eqnarray*}&&t=t/{T}_{o},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{T}_{o}={n}_{o}/({{gX}}_{o}).\end{eqnarray*}$$

Values used for the characteristic scales are X_o ∼ 10⁴ m, H_o ∼ 250 m, n_o ranging from 10³–10¹¹ Pa s (for comparison, n_H2O ∼ 10⁻³ Pa s), and gravity g = 3.7 m s⁻². In the process of collecting terms in the scaling of the governing equations, we ended up with two nondimensional parameters,

$$\begin{eqnarray*}&&e={H}_{o}/{X}_{o},\mathrm{and}\end{eqnarray*}$$

$$\begin{eqnarray*}&&b={{n}_{o}}^{2}/({{{gX}}_{o}}^{3}),\end{eqnarray*}$$

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) and (E) range from several hours to days for an S₈-dominated flow or many years for a halite-dominated 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⁻¹ on small slopes (Takagi 2010). The effective viscosity of various mafic lava flows on Earth is in the range of 10³–10⁶ Pa s (Chevrel et al. 2019; water has a viscosity of ∼10⁻³ Pa s at 20 °C). Ice glaciers deform by slow creep of dislocations and exhibit an effective viscosity of roughly 10¹² Pa s with typical flow rates of 0.1–1 m day⁻¹. The effective viscosity of ductile halite flow at ∼200 °C (Li & Urai 2016) is about 6 × 10¹⁵ Pa s. At 400 °C, halite viscosity drops to about 10¹⁴ Pa s.

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 millennia-spanning 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, volatile-enriched 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 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 finite-difference 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

$$\begin{eqnarray*}&&\displaystyle \frac{\partial T}{\partial t}={k}_{T}\displaystyle \frac{{\partial }^{2}T}{\partial {x}^{2}}\end{eqnarray*}$$

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⁻¹ K⁻¹, ρ = 3000 kg m⁻³, and c_p = 1000 J kg⁻¹ K⁻¹.

At the surface,

$$\begin{eqnarray*}&&-{K}_{T}\displaystyle \frac{\partial T}{\partial x}=(1-\alpha )S(t)-\sigma \varepsilon {T}^{4},\end{eqnarray*}$$

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⁻², P is the rotational period of 58 Earth days, ε is the emissivity - assumed to be 0.95, and σ is the Stefan–Boltzmann coefficient (5.6704 × 10⁻⁸ W m⁻² K⁻⁴). At depth (500 m here),

$$\begin{eqnarray*}&&-{K}_{T}\displaystyle \frac{\partial T}{\partial x}=G,\end{eqnarray*}$$

where G is the geothermal heat flux, assumed to be 20 mW m⁻². 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⁻¹ K⁻¹, but a lower value of 3 W m⁻¹ K⁻¹ 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⁻¹ (Figures 3 and 4; Breuer et al. 2007). If mantle thermal conductivity is as high as 4 W m⁻¹ K⁻¹, then the heat flux to the surface would be 0.020 W m⁻², or 20 mW m⁻², 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-of-magnitude 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 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, 2013, 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).