Buried Ice Deposits in Lunar Polar Cold Traps Were Disrupted by Ballistic Sedimentation

The NASA Artemis program will send humans to the lunar south polar region, in part to investigate the availability of water ice and other in situ resources. While trace amounts of ice have been detected at the surface of polar permanently shadowed regions (PSRs), recent studies suggest that large ice deposits could be stable below cold traps in the PSRs over geologic time. A recent study modeling the rate of ice delivery, ejecta deposition and ice loss from cold traps predicted that gigatons of ice could be buried below 100s of meters of crater ejecta and regolith. However, crater ejecta vigorously mix the target on impact through ballistic sedimentation, which may disrupt buried ice deposits. Here, we developed a thermal model to predict ice stability during ballistic sedimentation events. We then modeled cold trap ice and ejecta stratigraphy over geologic time using Monte Carlo methods. We found that ballistic sedimentation disrupted large ice deposits in most cases, dispersing them into smaller layers. Ice retention decreased in most cases, but varied significantly with the sequence of ejecta delivery, particularly from basin‐forming events. Over many model runs, we found that south polar craters Amundsen, Cabeus, and Cabeus B were most likely to retain large deposits of ice at depths up to 100 m, shallow enough to be detectable with ground‐penetrating radar. We discuss these findings in the context of the imminent human exploration activities at the lunar south pole.

The first model investigating ice and ejecta stratigraphy at the lunar poles was developed by Cannon et al. (2020).Through a Monte Carlo approach, they showed that ejecta blanketing of cold traps could preserve mining-scale ice deposits over geologic time.However, that work neglected the effects of ballistic sedimentation, the vigorous mixing of ejecta with local materials, which may rework and volatilize ice rather than strictly preserving it (Kring, 2020;Oberbeck, 1975;Weiss & Head, 2016).In this work, we seek to understand the effects of ballistic sedimentation on lunar polar ice and ejecta stratigraphy.In order to address this question, we developed a simple thermal model to account for volatilization in a particular ballistic sedimentation event.We then developed a Monte Carlo polar ice and ejecta stratigraphy model using the same framework of Cannon et al. (2020).In addition to ballistic sedimentation, we use the model to explore the effects of basin impacts, cometary impactors, and solar wind ice deposition.We predict ice retention at key cold trap locations in the Artemis exploration zone and discuss the potential for subsurface ice exploration near the lunar south pole.

Methods and Modules
Moon Polar Ice and Ejecta Stratigraphy (MoonPIES) is a Monte Carlo model designed to simulate lunar polar cold trap stratigraphy resulting from ice delivery, ejecta deposition, and ice removal at the lunar poles.Our model extends a previous model by Cannon et al. (2020) by introducing basin-scale impacts, ballistic sedimentation effects, latitude dependent ballistic hop efficiency, and cometary impactors.

Main Model
We developed the MoonPIES model to track ice and ejecta layering within permanent cold traps over lunar geologic time, recording ice delivered, ice lost, and ejecta deposited to target cold traps in 10 Myr intervals.We 10.1029/2022JE007567 3 of 23 limited our study to south polar cold traps found within large (D > 20 km) craters (Cannon et al., 2020).Each cold trap was modeled as a 1D column of unit area at the centroid of its cold trap area (Figure 5; Data Set S1 in Supporting Information S1).
At each time step following its formation age, a given cold trap stratigraphy column was updated in the following order: 1. Ejecta deposited from basins and craters (possible ballistic sedimentation) 2. Ice deposited (see processes, Figure 1) 3. Ice removed due to impact gardening (Section 2.2) We ran the full MoonPIES Monte Carlo model from 4.25 Ga to the present 10,000 times to generate a statistical distribution of possible cold trap stratigraphies.Parameters which were varied on each run included the ages of basins and polar craters (Figure 2) as well as the amount of impactor ice deposited in each time step (see Section 2.5).

Ice Loss by Impact Gardening
Ice loss from lunar cold traps is poorly constrained; however, Farrell et al. (2019) suggested that ice can be fully eroded from the upper 500 nm of a cold trap on ∼2,000 years timescales due primarily to the effects of micrometeorite bombardment.Furthermore, plasma sputtering and thermal effects may enhance this loss (Farrell et al., 2019).Using an analytical model of impact gardening, the constant churning of the regolith due to impacts, Costello et al. (2020Costello et al. ( , 2021) ) found that buried ice is exposed to near-surface loss processes more rapidly than it is buried.Therefore, impact gardening on long timescales (model time steps of 10 Myr ≫ 2 kyr surface stability), 4 of 23 provides ample opportunity for buried ice to be exhumed and subsequently vaporized, sputtered, or desorbed from the surface.
The previous ice stratigraphy model by Cannon et al. (2020) used impact gardening to estimate an ice loss rate of 10 cm per 10 Myr, the depth of the heavily mixed "in situ reworking zone" (Costello et al., 2018).In our model, we update this value to 9 cm per 10 Myr, with an increase prior to 3 Ga (Costello et al., 2021).Figure 6 of Costello et al. (2021) shows an excess reworking depth prior to 3 Ga of approximately 1% of the early impact flux relative to the present-day flux (Neukum et al., 2001).Therefore we increase the reworking depth by 0.9 mm per 10 Myr before 3 Ga, yielding 20 cm per 10 Myr in the first model time step (4.25 Ga).Since the in situ reworking zone is vigorously mixed to homogeneity (Costello et al., 2021), we assume all ice to this depth is exposed to surface loss at some point in the 10 Myr and therefore our model removes all ice within the upper reworking depth each time step.A consequence of this treatment is that any ejecta within the upper reworking depth preserves underlying ice according to its thickness (e.g., for a reworking depth of 9 cm, surface ejecta layers at least 9 cm thick will preserve underlying ice layers, while 5 cm of near-surface ejecta allows 4 cm of ice to be lost).In the pessimistic case where no ejecta preservation occurs, about 9 m of total ice thickness could be lost per Gyr in each of the last 3 Ga of model time, while up to 45 m of ice could be lost over the full 4.25 Gyr when accounting for the deeper reworking depth earlier than 3 Gyr (Figure S10 in Supporting Information S1).Our simple estimate of ice loss as a function of impact gardening does not account for partial mixing and redistribution of ice deeper than the in  Tye et al. (2015), Deutsch et al. (2020), andCannon et al. (2020).Lunar epochs were defined by the ages of the basins (Orgel et al., 2018).Asterisks indicate a basin that deposits ejecta into at least one of the modeled craters.

Ejecta Deposition
Ejecta deposition has two main functions in the model: destruction of buried ice through ballistic sedimentation (Section 2.4) and preservation of ice against impact gardening and surface loss processes (Section 2.2).We modeled ejecta emplacement similarly to Cannon et al. (2020) with the addition of basin events (Figure 2).We included 24 south polar complex craters (D > 20 km) and 27 basins (D > 300 km), each dated previously by crater counting methods (Data Sets S1 and S2 in Supporting Information S1; Deutsch et al., 2020;Tye et al., 2015;Orgel et al., 2018).We calculated ejecta thickness (t) as a function of distance from the crater center (r) and crater radius (R) using the scaling relationship from McGetchin et al. (1973), Kring (1995): 0.04 for simple craters, 0.14 for complex craters An exponent of −2.8 ± 0.5 was proposed by Fassett et al. (2011) for basin-size events, but that is within the uncertainty of Equation 1. Two caveats of using Equation 1to estimate ejecta thickness are as follows: (a) proximal ejecta may be overestimated due to the raised crater rim (Kring, 2007;Sharpton, 2014) and (b) distal ejecta is heterogeneously distributed (e.g., as discontinuous rays; Gault et al., 1974) due to instabilities in the ejecta curtain (Melosh, 1989).Because none of the craters or basins in our study overlap one another, no cold traps lie on or within 1 crater radius of another crater.Furthermore, we restricted ejecta deposition to threshold distances at which most impacts produce continuous ejecta (Melosh, 1989), setting the threshold distances to be 4 crater radii from the center of each crater, for consistency with Cannon et al. (2020), and 5 radii from the center of each basin (T.Liu et al., 2020;Xie et al., 2020).Basins that are near enough to deposit ejecta into a modeled south polar cold trap crater are indicated with asterisks in Figure 2.

Ballistic Sedimentation
Ballistic sedimentation was first formally discussed by Oberbeck (1975) and describes the process by which ejecta from a primary crater follows a ballistic trajectory and impacts the surface at high velocity, mixing with local material to form breccias.A side effect of ballistic sedimentation is heating of the mixed ejecta unit.Because we have not drilled into ballistic sedimentation breccias on the Moon, we use Earth analogs to constrain the effects of ballistic sedimentation, specifically, the Bunte Breccia Unit within Ries crater in Germany (Hörz et al., 1977(Hörz et al., , 1983;;Oberbeck, 1975).
Ballistic sedimentation only occurs when ejecta reach the surface with sufficient velocity to mix the target.
In the case of Meteor Crater (D = 1.25 km; see Kring, 2007, for a review), a simple crater on Earth, continuous ejecta is distributed to distances of about two crater radii beyond the crater rim.Material deposited at the outer edge of that ejecta blanket hit the surface with a velocity of about 11 m/s, which caused some radially outward skating across the landscape but no significant erosion and mixing with substrate materials or heating.Conversely, the larger Ries crater (D = 24 km) has two distinct ejecta units.It contains a polymict breccia with fragments of solidified impact melt (known as suevite) that has components shocked to >50 GPa and depositional temperatures between 500 and 900°C (Kring, 2005, and references therein).The underlying Bunte breccia unit represents the bulk of the ejecta and has components shocked to <10 GPa and was deposited with no significant increase in temperature.The Bunte Breccia is the ballistically emplaced unit at the Ries Crater (Hörz et al., 1977(Hörz et al., , 1983)).The Bunte Breccia extended from near the crater rim to at least 36 km from the crater center, about three times the basin radius (see Kring, 2005, for a review).We therefore expect the onset of ballistic sedimentation to occur at an ejecta kinetic energy intermediate to those at Meteor and Ries craters.
Since ballistic sedimentation is capable of introducing significant heating and mixing to target materials, we modeled it as an ice loss process as a function of ejecta kinetic energy, temperature, and depth of mixing.

Onset of Ballistic Sedimentation
To estimate a threshold kinetic energy for the onset of ballistic sedimentation, we computed the relevant ejecta energy of the Bunte Breccia deposits at the Ries impact structure (Kring, 2005).We computed ejecta mass per unit area as the product of ejecta density and thickness, assuming an ejecta density of ρ = 2,700 kg/m 3 corresponding to the locally ejected Malmian limestone (Bohnsack et al., 2020) and calculating thickness using Equation 1 (R Ries = 12 km).Ejecta velocity at impact was computed using the ballistic formula for a spherical body (Vickery, 1986), derived from the half-angular distance of travel, ϕ = r/2R p , which is related to velocity (v) and the ejection angle (θ) by the following equation: where g is gravitational acceleration and R p is planet radius.Solving for velocity gives the following equation: (3) Ballistic kinetic energy is then given by KE = mv 2 and is a function of r, R p , and θ.Using the 4 crater radii extent of our ejecta deposits and Equation 3 for the Earth (g = 9.81 m/s 2 ; R p = 6,371 km) and assuming the most likely ejecta angle of θ = 45° (Shoemaker, 1962), we find that the minimum kinetic energy that produced ballistic sedimentation at the Ries impact structure was ∼1,500 MJ/m 2 (Figure 3c).
For comparison, we repeated the above calculation for the 1 km Meteor Crater where ballistic sedimentation has not been observed (see Kring, 2007, for a review).The maximum ejecta kinetic energy was ∼10 MJ/m 2 at Meteor crater, indicating the onset of ballistic sedimentation is between 10 and 1,500 MJ/m 2 .
The lunar case for ballistic sedimentation as a function of crater size is made by Oberbeck (1975) who observed the onset of hummocky textures in the continuous ejecta of craters larger than 4 km in diameter.Repeating the kinetic energy calculation for the Moon assuming an anorthositic target (ρ = 2,700 kg/m 3 ; g = 1.624 m/s 2 ; R p = 1,737 km), we found the kinetic energy for Meteor Crater was smaller than the 4 km lunar crater (Figure 3).If the hummocky textures observed by Oberbeck (1975) are indicative of ballistic sedimentation, we would expect the Meteor Crater impact to be energetic enough to produce these deposits.However, since the smallest known crater to produce ballistic sedimentation deposits is the Ries impact, we use it as our conservative model threshold.The derived ∼1,500 MJ/m 2 at 4 crater radii corresponds to lunar craters D ≥ 20 km, and therefore we model ballistic sedimentation for impacts larger than 20 km which encountered a modeled south polar cold trap crater (Figure 2) within 4 primary crater radii.Although ballistic sedimentation events by craters 4-20 km are not captured in our model, those events would primarily influence shallower and younger ice deposits which are not the primary focus of this work (see Section 4.2).

Ballistic Sedimentation Depth
We modeled the depth of influence of ballistic sedimentation events using the local mixing model introduced by Oberbeck (1975) and updated by Petro and Pieters (2004).The local mixing ratio, μ, of local material to ejected material was modeled as a function of distance traveled, r in km (Oberbeck, 1975), with an adjustment for μ > 5 (Petro & Pieters, 2006): We then defined the ballistic sedimentation depth (δ) as the product of ejecta thickness (Equation 1) and mixing ratio (Petro & Pieters, 2004): The mixing ratio can also be expressed as the fraction of ejecta relative to total (ejecta and local) material as follows: We parameterize ice loss to ballistic sedimentation as a function of f ejecta and the temperature of the incoming ejecta.

Fraction of Ice Volatilized
We used a 1D heat flow model to derive the fraction of local material volatilized (heated beyond a constant ice stability temperature), given the amount of ejecta delivered (f ejecta ) and the initial ejecta and target temperatures.
The 1D heat flow model assumes vigorous mixing and rapid equilibration, such that the ejecta and local material primarily exchange heat through conduction rather than radiation (Carslaw & Jaeger, 1959;Onorato et al., 1978).
We then solve the 1D heat flow equation using the forward Euler method (Euler, 1792) and track the maximum temperatures encountered during equilibration.We summarize the model in Equations 7 and 8: where K is the thermal diffusivity, κ(T) is thermal conductivity (see Equation A4 in Hayne et al., 2017), ρ = 1,800 kg m −3 is regolith density, C p (T) is the heat capacity (see Table A1 in Hayne et al., 2017), Δt = 1 ms is the time step, Δχ = 10 μm is the spatial scale, and i, n are the spatial and time steps, respectively.
To compute the fraction of water ice volatilized ("volatilized fraction" hereafter), we chose a surface ice stability temperature of 110 K (Fisher et al., 2017;Hayne et al., 2015).A regolith overburden may increase the ice stability temperature (Schorghofer & Aharonson, 2014); however, we do not track the depths of particles in the model and using the surface ice stability provides an upper bound on ice mobilized during the rapid ballistic sedimentation mixing process.We initialized the 100 element zero-dimensional model (closed 1D loop tracking only heat exchanged conductively between participating grains) to a nominal surface PSR temperature of 45 K (typical surface temperatures of lunar polar cold traps; Paige et al., 2010) and randomly assigned a chosen number of ejecta elements (dictated by f ejecta ) a particular initial ejecta temperature.We varied f ejecta and ejecta temperature from 0% to 100% and 110 to 500 K, respectively.We ran the thermal model until equilibration (i.e., all elements within 1 K of each other), or until all elements exceeded 110 K (volatilized fraction = 1).We defined the volatilized fraction as the fraction of local elements that exceeded 110 K in any time step.We randomized the initial ejecta positions and ran the model 50 times for each f ejecta and ejecta temperature combination, reporting the mean and standard deviation volatilized fractions of all runs (Figure 4; Data Sets S3 and S4 in Supporting Information S1).

Ice Lost Due To Ballistic Sedimentation
To predict ice loss in a particular ballistic sedimentation event, we estimated the ejecta and target cold trap temperatures at the time of deposition.Hydrocode simulations by Fernandes and Artemieva (2012) indicated that ejecta temperatures increase with distance from a basin impact due to shock heating, but primarily beyond the 4-5 radius distances modeled here.Proximal ejecta temperatures were much more sensitive to the choice of subsurface thermal profile from about 260 K in the "present-cold-Moon" scenario to 420 K for the "ancient-hot-Moon" (Artemieva & Shuvalov, 2008).For a conservative treatment, we chose 260 K as the ejecta temperature of basin impact ejecta (basin D ≥ 300 km), but note that there is little change in our predicted volatilized fraction from 260 to 420 K (Figure 4).For smaller polar complex craters (20 km ≤ D ≤ 110 km), we chose an ejecta temperature of 140 K, a typical subsurface polar regolith temperature (Feng & Siegler, 2021;Vasavada et al., 1999).
For each ballistic sedimentation event, we then retrieved f ejecta (Equation 6) and ejecta temperature (Figure 4) and computed the ballistic sedimentation depth, δ (Equation 5).We assumed that the derived volatilized fraction of ice was volatilized in each layer within δ of the surface.All volatilized ice is assumed to be lost from the stratigraphy column in the ballistic sedimentation event since each column tracks only the ingress and egress of ice to a permanently shadowed macro cold trap.We discuss the implications of this simple treatment of ballistic sedimentation ice removal in Section 4.2.If multiple ballistic sedimentation events occurred in a single time step, 10.1029/2022JE007567 9 of 23 they were applied in δ order from smallest to largest.The ballistic sedimentation events in our model produced δ values ranging from meters to multiple kilometers depending on the size of primary impact and distance to a particular cold trap (Figure 3).Our first order approximation of ballistic sedimentation effects allows us to assess incoming ejecta as a removal process; however a method which more precisely mixes and redistributes ice may be warranted in future work.

Impact Ice Delivery
To model ice delivery to the lunar poles by impacts, we first divide all possible impactors into six size regimes and two classes: hydrated asteroid and comet.Hydrated asteroids in Regimes B-E were modeled consistently with Cannon et al. (2020) as 24% hydrated C-types with 10% water content by mass, and adopting consistent fluxes, size-frequency distributions, and crater scaling laws, as summarized in Table 1 and Figure S10 in Supporting Information S1 (Brown et al., 2002;Grün et al., 2011;Neukum et al., 2001;Ong et al., 2010).Additionally, we introduce basin impactor (the new Regime F) and cometary impactor (across all regimes) contributions to polar ice for the first time.Note.Regimes A-E are defined following (Cannon et al., 2020).Regime F, representing basin impactors, was added for this work.a MoonPIES only.

Micrometeorites and Comets
We modeled all micrometeorites (Regime A) as cometary based on dynamical models that suggest comets dominate the smallest impactors at 1 AU (Oberst et al., 2012;Pokorný et al., 2019;Suggs et al., 2014).For larger impacts, we assumed 5% are cometary, which is a conservative estimate within a typical range of 5%-17% (Joy et al., 2012;J. Liu et al., 2015).To predict the total cometary mass delivered, we used the same scaling laws of Regimes B-F (Table 1).
We assigned each comet a random velocity from a bimodal velocity distribution (Figure S2 in Supporting Information S1) and assumed each comet is 20 wt% water, which has been measured for comet 67P/Churyumov-Gerasimenko (Fulle et al., 2017), but is on the low end of historical estimates, for example, 50 wt% (Whipple, 1950), with a conservative density of 600 kg m −3 (the density of Comet Shoemaker-Levy 9; Asphaug & Benz, 1994).
We model the micrometeoritic flux rate at 10 6 kg/yr (Grün et al., 2011), scaled by the lunar chronology function, and assuming complete vaporization and consistent cometary hydration of 20 wt%.We note that the micrometeorite flux rate and hydration depend on the sizes and source populations of parent comets which have evolved through time (Pokorný et al., 2019;Snodgrass et al., 2011).Although we do not model comet populations through time, we use conservative estimates of comet parameters for this work and show the relative insensitivity of the model to very icy comets in a 100 wt% comet hydration scenario (Figure S6 in Supporting Information S1).

Basins
Basins included in our model (Figure 2; Data Set S2 in Supporting Information S1) were designated as hydrated C-type, cometary, or neither at the same probability as other impactors.If a modeled basin impactor was icy, it delivered ice in the time step nearest to its age, randomized within model age uncertainties for each run.Although absolute ages and age uncertainties are debated for several basins, we drew all basin ages from the same source, Orgel et al. (2018), noting that the sequence of basin ages is more important for generating an accurate stratigraphy than precise absolute ages.We scaled each basin main ring diameter to its transient diameter using the scaling laws from Croft (1985) and then to impactor diameter following Johnson et al. (2016).Volatile content and retention rates were predicted consistently with other hydrated asteroids and comets in the model.

Volatile Retention
The fraction of ice retained by the Moon in a hypervelocity impact is primarily a function of the impact velocity (Ong et al., 2010).We derived a simple power law fit to retention rates computed in these impact simulations (1.66 × 10 4 v −4.16 , where v is velocity in m/s; Figure S3 in Supporting Information S1).For consistency with Cannon et al. (2020), we retained the assumption that <10 km/s impacts result in 50% volatile retention for asteroid impactors.

Ballistic Hopping
At each time step, global ice delivery was converted to a local cold trap ice mass by employing a "ballistic hop efficiency," defined as the fraction of global ice that comes to rest in a particular cold trap via ballistic hop random walks.We took the fraction of total water in the ballistic hop simulations by Moores (2016) and normalized by cold trap area.For cold traps not modeled in Moores (2016), we made a conservative estimate recognizing that ballistic hop efficiency is related to latitude.For Slater, we took the average of the nearest latitude craters, Shoemaker, de Gerlache, and Sverdrup.For craters north of Faustini (87.2°S), we set the ballistic hop efficiency to that of Faustini, recognizing that if the latitude trend holds then this would underestimate the amount of ice transported to and therefore retained within these cold traps in the model (Figure S1 in Supporting Information S1; Data Set S5 in Supporting Information S1).We note that recent work has called into question the present-day lateral transport of ice (Hodges & Farrell, 2022), while other studies suggest impact-induced transient atmospheres may also play a role in ice delivery (Nemtchinov et al., 2002;Stewart et al., 2011).While we use the transient atmosphere formulation for volcanic ice delivery (Section 2.6), we retain the ballistic hop method for impact ice delivery for ease of comparison with Cannon et al. (2020).

Volcanic Ice Delivery
We modeled volcanic ice delivery via a transient atmosphere that deposits ice in polar cold traps (Aleinov et al., 2019).Volatiles are deposited at a rate predicted by transient atmosphere simulations by Wilcoski et al. (2021), who found 26% of erupted H 2 O is able to be deposited in south polar cold traps when accounting for atmospheric escape and sublimation.This treatment of deposition ignores any effect of ballistic hopping as the deposited H 2 O from a transient atmosphere can only persist on the poles and is very quickly sublimated away elsewhere on the surface (Wilcoski et al., 2021).We used model estimates of Needham and Kring (2017) for total H 2 O outgassed from mare volcanic provinces over time.We converted volatile H 2 O to ice deposited in the style of Cannon et al. (2020).Although Head et al. (2020) presented smaller estimates of outgassed volatile mass, our model is insensitive to this choice as neither deposits more ice than is removed by impact gardening in a given time step (Figure 6, Figure S10 in Supporting Information S1).

Solar Wind H + Deposition
We included a treatment of solar wind H + as a possible source of water in polar cold traps (Arnold, 1979).We used solar wind-derived H 2 O mass flux of 2 g/s H 2 O (Benna et al., 2019;Housley et al., 1973).We note this may be an overestimate since Lucey et al. (2020) predicted about 1/1,000 of the 30 g/s H + suggested by Hurley et al. (2017) would be converted to H 2 O. Additionally, historical solar wind H + flux may have been lower when the sun was fainter (Bahcall et al., 2001).However, less than 1 cm of ice is expected per model time step given the 2 g/s H 2 O rate, which is less than is removed by impact gardening each time step, making our model insensi tive to a more precise treatment of solar wind ice delivery (Figure 6, Figure S10 in Supporting Information S1).

Randomness and Reproducibility
To simulate the delivery of ice and ejecta to target cold traps, MoonPIES uses Monte Carlo methods to vary the timing and abundance of impacts through lunar geologic history.Consistent with the previous model by Cannon et al. (2020), ice delivery was driven by impactor size and impactors forming craters smaller than 1.5 km diameter were treated as bulk populations, with fluxes drawn from the literature (see Table 1).However, larger impactors were modeled individually with compositions (hydrated asteroid, unhydrated asteroid, or comet) and speeds each randomly drawn, affecting the ice delivered and retained in a given time step (see Section 2.5).Asteroid speeds were drawn from a Gaussian distribution (μ = 20 km/s, σ = 6 km/s).Comet speeds were drawn from a bimodal distribution of two Gaussians to simulate the Jupiter family and long period comet populations (μ JFC = 20 km/s, σ JFC = 5 km/s; μ LPC = 54 km/s, σ LPC = 5 km/s) (Figure S2 in Supporting Information S1; Chyba, 1991;Jeffers et al., 2001;Ong et al., 2010).We assume that Jupiter family comets are seven times as likely as long period comets, though the precise ratio and its evolution over lunar geologic time are not well constrained (Carrillo-Sánchez et al., 2016;Pokorný et al., 2019).
To model individual crater and basin events that deliver ejecta to cold traps, we randomly drew ages from the published crater count model ages (Figure 2).Ages were drawn from a truncated Gaussian centered on the absolute model age with standard deviation being half of the model age uncertainty (or the average of the upper and Figure 6.Mean ice deposited in each time step per unit area across the south polar region for each lunar geologic era (log scale, averaged over 100 runs).Whiskers denote the maximum and minimum ice delivered in a particular time step.Bars show total ice delivered from all sources (blue) as well as ice originating from nonbasin impactors (orange), basin impactors (green), volcanic outgassing (red), and solar wind (yellow).Shaded regions represent upper bounds on ice loss due to ballistic sedimentation (pink) and impact gardening (gray).In practice, ice loss will depend on the timing of stochastic ballistic sedimentation events and stratigraphy of each cold trap at each time step, since surface ejecta layers may partially or fully preserve underlying ice from loss (to visualize stochasticity and cumulative deposition, see Figure S10 in Supporting Information S1).
lower bounds, if asymmetric).The minimum and maximum ages are fixed at the upper and lower age uncertainty.Ice delivery was not modeled for individual south polar craters since their contributions are already included in the complex crater population.However, basins are assigned a probability of being asteroidal or cometary and a random impact speed to determine ice delivery in a consistent manner with other impactors.Ejecta deposition and basin ice delivery then occurs during the nearest time step in the model to the assigned ages for these individual events.
For reproducibility, MoonPIES generates a configuration file specifying all parameters used to run the model as well as the random seed, allowing a given model run to be reproducible.Monte Carlo results presented here were run on MoonPIES v1.1.0with random seeds 1-10000 (see Tai Udovicic et al., 2022a, for full documentation).

Results
The ice layering trends predicted by the MoonPIES model were broadly consistent with Cannon et al. (2020) when ballistic sedimentation was excluded from the model (e.g., the oldest cold traps often retained "gigaton" ice deposits, Section 3.2).However, when accounting for ballistic sedimentation, such gigaton ice deposits were disrupted and the quantity of ice retained was reduced.We also found a location dependence of ice retention, with cold traps nearest to the south pole (latitude <−88°) retaining less ice in general when compared to cold traps at greater latitudes (Figure 9).We ran the model 10,000 times to generate a distribution of ice retention in each cold trap, with the results of the Monte Carlo approach shown in Figures 7 and 9. Figure 7 compares the stratigraphic columns generated for Faustini, Haworth, and Cabeus craters over three different Monte Carlo model runs, while Figure 9 shows the distribution of total ice thickness in each of the 10 modeled stratigraphic columns across all 10,000 model runs (Data Set S6 in Supporting Information S1).

Effects of Ballistic Sedimentation
Ballistic sedimentation reduced the amount of ice retained within most cold traps modeled in this study.Figure 8 depicts a single run with and without the effects of ballistic sedimentation, while Figure 9 depicts total ice thickness retained for 10 cold traps across 10,000 model runs with and without the effects of ballistic sedimentation.Cold traps were grouped by similar age and sorted by latitude within each group.
When ballistic sedimentation is considered, ejecta from nearby craters and impact basins effectively removed preexisting ice within cold traps, resulting in less preserved average ice thickness across all cold traps.Ballistic sedimentation did not remove ice entirely but rather reduced the thickness of ice deposits at and near the surface at the time of ejecta implantation.Although impacts that form basins and complex craters both have the potential to remove ice due to ballistic sedimentation, basin-sized impacts tended to be more effective at melting ice through ballistic sedimentation than complex polar craters.Complex polar craters tended to produce ejecta lacking sufficient kinetic energies at the distances required to reach nearby polar cold traps, effectively caused a net preservation effect rather than a net ice removal effect.The influence of ballistic sedimentation was most notable for the oldest cold traps in our sample, Haworth, Shoemaker, and Faustini.The oldest cold traps retained large ice deposits in nearly all runs without ballistic sedimentation.After ballistic sedimentation effects were included, ice retention was reduced and more variable, with cold traps retaining about a tenth of former median total ice thicknesses (Figure 9).
Similarly, Nectarian and Imbrian cold traps saw overall declines in ice retention when including the effects of ballistic sedimentation, though not as dramatically as observed in the Pre-Nectarian cold traps, which follows from the decline in basin-forming impacts during these eras.Eratosthenian cold traps formed after the majority of ice was delivered and most large ballistic sedimentation events occurred and therefore retain little ice regardless of ballistic sedimentation (Figure 9).

Gigaton Ice Deposits
We report large ice deposits retained in the model following the "gigaton" terminology, referring to ice deposits that would exceed 10 9 metric tons of ice if they filled a cold trap with a surface area of at least 100 km 2 (Cannon et al., 2020).Assuming a density of ice of ∼1,000 kg/m 3 , such deposits must be larger than 10 m thick.Pristine single-layer deposits exceeding 10 m were rare in most cold traps and absent from others except when early large ice delivery events occurred in a particular run (Figure 8).We investigated the likelihood of retaining layers of a given thickness at depth over the 10,000 model runs (Figure 10).When excluding effects of ballistic sedimentation, gigaton layers (>10 m) emerged at depths of about 100 m.By contrast, when ballistic sedimentation was implemented, layers rarely exceeded 10 m and were most commonly <1 m thick (Figure 10).
We also assessed the total ice retention by each cold trap at all depths (Figure 9).Without ballistic sedimentation, we found that the median total ice exceeded 10 m "gigaton" thickness for all Nectarian and older craters.However, when we accounted for ballistic sedimentation, all craters declined in total ice retention, with medians near or below 10 m.For Shoemaker, Idel'son L and Amundsen, this decline marked a shift from >75% of runs producing gigaton levels of ice to <25% exceeding the 10 m threshold.However, other Nectarian and older cold traps retained 10 m of ice in about 50% of runs.The youngest cold traps retained the least ice: Imbrian cold traps Slater and Sverdrup retained at least 1 m of ice in about 15% and 10% of runs, respectively, while Eratosthenian cold trap Wiechert J exceeded 1 m in only 0.5% of runs.Shackleton never retained more than 1 m of ice.In summary, while single layer gigaton deposits were rare, volumetric gigaton deposits occurred for most older cold traps in about 50% of runs.

Gardened Layer
The gardened layer (referred to as gardened mantle deposits in Cannon et al. ( 2020)) that we observed near the surface of most stratigraphy columns ranged in size depending on cold trap location and age.We defined the gardened layer as any portion of the column influenced by ballistic sedimentation or impact gardening and containing less than the 100 m "gigaton" ice thickness threshold defined above.For any given cold trap, the total thickness of the gardened regolith zone varied significantly across Monte Carlo runs. Figure 9 demonstrates that the distribution of total ice thickness for 10,000 model runs varies from 1 to 100 m for Pre-Nectarian and Nectarian cold traps, while younger cold traps tend to retain between 1 and 10 m.S7 in Supporting Information S1).We note that the cold trap ages and stratigraphic sequence of a particular run differs due to random variation in ejecta ages and ice delivery.Runs A and B show typical columns for all three cold traps.In run C, the large ice layer in Faustini and Haworth is the result of a basin-scale cometary impact, while the absence of layering in Cabeus indicates that it formed later than all possible ejecta sources. 10.1029/2022JE007567 14 of 23 We also observed a slight trend with latitude, finding that median total ice in the gardened layer declines toward the South Pole.In particular, Faustini, Cabeus B, and Cabeus retained the thickest gardened layers with median total ice of ∼20 m.The smallest gardened layers occurred in the youngest and most poleward cold traps, namely Shoemaker, de Gerlache, Slater, Sverdrup, Wiechert J and Shackleton, each retaining <10 m median ice thickness.Wiechert J and Shackleton were the only cold traps which retained <1 m of total ice in the majority of runs and exceeded this threshold in only 10% and 1% of runs, respectively.
We contrasted our results with 10,000 model runs which did not account for ballistic sedimentation (Figure 9).As expected, the total amount of ice retained in each cold trap was greater when ballistic sedimentation was excluded.In addition, thick, pristine ice layers were disrupted by ballistic sedimentation, causing the gardened regolith to be much more extensive in this work relative to the previous model by Cannon et al. (2020).S8 in Supporting Information S1).When ballistic sedimentation was accounted for, large pure ice layers were lost from the base of the oldest cold traps (Faustini, Haworth, and Shoemaker).Shallower layers retain a similar ice % in both cases.Although this is a representative outcome, it should be noted that the absolute quantity of ice and stratigraphic sequence changes from run to run (Figures 7-9).
In addition to total ice, we also illustrate the possible distributions of ice with depth for Faustini, Haworth, Amundsen, Cabeus, de Gerlache, and Slater craters in Figure 10 using a kernel density estimation (KDE) contour plot of ice thickness versus depth.The contours represent the number of times an ice layer of the corresponding thickness was present at the corresponding depth over 10,000 model runs, both with (blue) and without (brown) ballistic sedimentation effects.Gray boxes in the top two plots indicate the gigaton deposit zone.In this figure, we have indicated two possible excavation depths from the LCROSS impact for Cabeus (Luchsinger et al., 2021;Schultz et al., 2010).
For Faustini, Haworth, and Amundsen, the three oldest cold traps, the blue ballistic sedimentation contours are shifted to the left relative to the brown contours, which represent model runs with no ballistic sedimentation (Figure 10).The shift in the blue contours represents ice deposits that have been disrupted and reduced in thickness.Individual ice layers in the gigaton deposit zone, indicated by gray boxes, only occur in some model runs without ballistic sedimentation, and only for Faustini and Haworth craters.The contour lines for Cabeus, de Gerlache, and Slater are less affected by ballistic sedimentation.Faustini and de Gerlache both contain ice deposits within the uppermost 6 m of regolith in some model runs, potentially making these craters high priority for missions with depth sensitivity greater than 1 m.However, the depth and temporal resolution of the MoonPIES model limits its ability to predict surface expression of ice.
In Figure 11, we present a boxplot of ice concentration for all cold traps within the uppermost 6 and 10 m below the surface.The boxes denote the first and third quartiles, while the whiskers denote the 99th percentile, and the individual points represent outliers above the 99th percentile.At least 75% of model runs predicted no ice retention in the uppermost 6 m for all cold traps except for Faustini, Amundsen, and de Gerlache.For comparison, the ∼5% concentration measured during the LCROSS impact into Cabeus crater (Colaprete et al., 2010) is indicated with a dashed line in Figure 11.None of the cold traps met or exceeded this threshold in the upper 6 m in >75% of model runs (Faustini and de Gerlache craters were most likely at 13% and 6% of model runs, respectively).Therefore, the MoonPIES model may underpredict ice near the surface, which we discuss in Section 4.2.However, we observe greater total ice retention at depths of 100 m, with Cabues B, Amundsen, and Cabeus cold traps exceeding 5% ice in the upper 100 m in at least 25% of model runs.We discuss the implications of these results in the following sections.

Gigaton Ice Deposit Distribution
Previous work by Cannon et al. (2020) described deposits containing tens to hundreds of meters of ice as "gigaton" deposits.When accounting for ballistic sedimentation in the MoonPIES model, we found that single-layer gigaton deposits no longer form (Figure 10).Instead, modeled stratigraphy columns are more likely to contain thinner layers that have been disrupted and reduced in size (Figure 10).It was also rare for the total ice thickness retained in our model to exceed 10 m throughout the stratigraphy column (Figure 9).This typically occurred when large volumes of ice were delivered by large icy impactors.Ice retention also depended on the precise sequencing of cold trap formation and ballistic sedimentation events.

Ice Disruption in Gardened Layers
Our model predicted that most ice layers stored within polar cold traps would have been disrupted by either impact gardening and/or ballistic sedimentation.While impact gardening affects the upper centimeters across all of model time, the effects of ballistic sedimentation are primarily localized to the Imbrian and earlier time periods and disrupts up to 10s of meters (Figure 6).We found that cold traps farther from the pole were more likely to retain thicker ice deposits, due to the MoonPIES model treatment of ballistic transport to the south polar region, wherein the majority of ballistically transported ice comes to rest in the first cold trap encountered where layer thickness exceeds 10 m.For Cabeus, horizontal lines denote LCROSS maximum excavation depth (Luchsinger et al., 2021;Schultz et al., 2010).(Moores, 2016).We note that our model only tracks deposition of ice by primary delivery processes and does not track the redeposition of ice that is lost from a particular cold trap, nor the interchange of ice between our target cold traps and seasonal or micro cold traps (Hayne et al., 2021;Kloos et al., 2019).In addition, ice volatilized during a ballistic sedimentation event may not fully escape the stratigraphy column.We tested the sensitivity of our model to ice lost during ballistic sedimentation events and find that at 50% ice lost, ballistic sedimentation still degrades the oldest, largest ice deposits, but becomes less effective for Nectarian and younger cold traps (Figure S5 in Supporting Information S1).Therefore, the gardened layers are conservative estimates of total ice storage.Secondary ice mobility processes such as redeposition and thermal pumping (Schorghofer & Aharonson, 2014;Schorghofer & Williams, 2020) could result in larger ice deposits than those modeled here, particularly in the case of near-surface expression of ice deposits.
In the MoonPIES model, the expression of near-surface ice in a gardened regolith was primarily determined by the thickness of the final ejecta deposition event(s) or large quantities of ice delivered by impacts near the end of a model run.The lack of surface ice observed follows from our treatment of impact gardening, which typically removed more ice than the average deposition from all sources in any given time step post-Nectarian (Figure 6, Figure S10 in Supporting Information S1).It is worth noting that reworking depth is a nonlinear function of time (Costello et al., 2021) and therefore ice loss due to impact gardening is sensitive to the chosen model time step.Shorter time steps increase the total number of gardening events, increasing the total ice loss in the pessimistic case (no ejecta preservation), albeit from shallower and easier to preserve reworking depths.For this reason, our model is more predictive of deeper ice layers rather than those directly at the surface.The surface expression of ice is therefore unconstrained by the MoonPIES model, as seen in Figure 11.MoonPIES also does not consider redeposition of lost ice into neighboring cold traps, thermal pumping, or other ice modification processes that may affect surface expression of ice.
Ice was most commonly retained at the surface in our model was when large icy impactors delivered ice very recently.Recent ice delivery has been invoked to explain the abundance of polar ice detected at Mercury and Ceres (Chabot et al., 2018;Moses et al., 1999;Platz et al., 2016), but not at the Moon (Campbell et al., 2006;Li et al., 2018;Neish et al., 2011).Instead, surface ice is more likely the result of present-day production of ice by the solar wind on timescales shorter than 10 Myr (Benna et al., 2019) or surface-subsurface exchange (e.g., through thermal pumping Schorghofer & Aharonson, 2014) not modeled here.Our model is therefore most relevant for subsurface ice exploration (e.g., depths of 3 m and greater).

The Nature of Buried Ice Below Polar Cold Traps
Our model aggregates ice delivery to the south pole from solar wind, volcanic outgassing, and impact delivery.
Using published rates of solar wind deposition (Benna et al., 2019;Housley et al., 1973) and volcanic outgassing (Needham & Kring, 2017), the ice deposition rate in all cold traps is many orders of magnitude smaller than the impact gardening rate at any point in time (Figure 6, Figure S10 in Supporting Information S1).Therefore, our model predicts that large buried ice deposits would be primarily sourced from impacts, consistent with the conclusions of Cannon et al. (2020).
Impactors in this work included contributions from both asteroids and comets, an update from previous work considering only asteroids (Cannon et al., 2020).We retained the asteroid ice delivery assumptions of Cannon et al. (2020), that is, 24% of asteroids were hydrated C-types with velocity-dependent volatile retention (Ong et al., 2010), except for velocities ≤10 km/s in which case heating only volatilized half of the carbonaceous projectile (Figure S3 in Supporting Information S1; Svetsov & Shuvalov, 2015).We assumed a constant comet proportion of 5% as a conservative estimate of typically reported values (5%-17%; Joy et al., 2012;J. Liu et al., 2015).At this proportion, comet ice contributions were comparable to asteroid ice due to the greater ice concentration, despite lower retention on average due to greater impact speeds (Ong et al., 2010).We also allowed basin impactors to be icy at the same rates as other impactors.On average, basins contributed about 2 m of ice per basin-era time step (1 m each from asteroidal and cometary basins; Figure S4 in Supporting Information S1).However, as stochastic events, many runs delivered no basin ice, while rare large basin impactors in other runs exceeded all other ice sources.In addition, basin ice delivery events produced deep early ice layers which were more likely to be retained over geologic time.Basin ice abundance and retention depended heavily on the sequence of basin events, their impactor composition and randomly assigned speed.Improved constraints on the fraction of cometary and asteroidal impactors over lunar history, their ice delivery mechanisms, as well as the ages and composition of basin impactors, would dramatically improve our understanding of deep water ice deposits at the lunar poles.

Implications for Lunar Ice Exploration
The MoonPIES model explores ice delivery, retention, and removal over geologic timescales.However, human exploration occurs on human timescales, during which short term ice deposition and removal can occur.These short term ice behaviors can lead to surface expressions of ice that are not captured by our model.Additionally, the true ice distribution is the result of a long stochastic history that we can only partially constrain, as illustrated by the variance in possible outcomes over 10,000 model runs (Figure 9).In particular, the bimodal distributions caused by uncertainty in the precise sequencing in crater and basin formation events indicate that precise sequencing of crater and basin ages is critical to the ability to precisely model and predict the thickness of ancient ice deposits.
Our model predicts that the disruption of ice by impact gardening and ballistic sedimentation would cause large coherent ice deposits at depth to be unlikely.We found that the gigaton deposits observed in Cannon et al. ( 2020) would be rare and, if present, likely to be disturbed and present only in incoherent layers.Therefore, well-defined contacts between lithic and icy layers are unlikely to be detected with radar, consistent with inconclusive space-borne radar studies (Campbell et al., 2006;Patterson et al., 2017;Spudis et al., 1998).If the LCROSS impactor excavated 10 m of material, it could have sampled a region that MoonPIES predicts could be populated by ice; however, if it only excavated 6 m of material, it would have sampled only the surface expression ice.Deposits beyond 6 m may have formed a layer of harder material that prevented excavation, as suggested by Luchsinger et al. (2021).
Additionally, ground-penetrating radar (GPR) may present an opportunity to probe for ice layers beneath cold traps with deeper penetration depths than orbital radar (Heggy et al., 2011;Kring, 2007Kring, , 2020;;Nunes et al., 2018;Ohtake et al., 2021;Richardson et al., 2020;Shoemaker et al., 2022;Sowers et al., 2022).GPR may also allow thinner ice layers to be detected by using higher frequencies than orbital radar.Although our predictions indicate that fully coherent thick ice layers are rare, changes in dielectric properties or partially preserved layers of ice may be observable with GPR.While Faustini, Haworth, and Shoemaker retained similar quantities of ice as other cold traps over all depths (Figure 9), most ice was concentrated near the base of their columns (Figure 10).The most valuable targets for radar assuming penetration depth of at least 100 m (Fa, 2013;Yuan et al., 2021) would be Amundsen, Cabeus, and Cabeus B, which each retained >3 m of ice in half of model runs and >5 m of ice in 25% of model runs (Figure 11).Faustini, Haworth, and Idel'son L, which each retained ≥1 m of ice in half of model runs may also be targets of interest.
The Artemis exploration zone is centered on the lunar south pole, which lies on the rim of Shackleton crater.Our model does not predict ice retention in Shackleton crater in the vast majority of model runs due to its relatively recent formation age and proximity to the south pole, consistent with orbital observations (Haruyama et al., 2008;Zuber et al., 2012).However, the absence of large subsurface ice layers in our model does not preclude the discovery of ice near the surface of Shackleton crater.Ice redistribution or a recent icy impactor could result in near-surface ice in the Shackleton cold trap.If a significant quantity of ice was discovered at depth below Shackleton, it would suggest that our model underestimates subsurface ice storage and could indicate that other cold traps may also store ice more efficiently than predicted here.Future exploration of south polar cold traps would therefore provide crucial constraints on our understanding of recent and historical ice delivery, as well as the potential for geologic deposits of ice at depth.

Conclusion
Understanding the location, quantity, and form of buried ice is critical for future mission planning.Volatile-bearing impacts are thought to be the main source of polar ice, while ejecta from impact craters may preserve ice deposits over geologic time.However, impact crater ejecta could mix and volatilize ice through ballistic sedimentation.
We developed a thermal model to predict ice loss due to ballistic sedimentation.We applied our findings to a Monte Carlo polar ice and stratigraphy model and determined that ballistic sedimentation disrupts "gigaton" style deposits reported by Cannon et al. (2020).Ice deposits in our model had smaller volume and layer sizes, particularly for older and deeper modeled ice layers.
We applied our model to cold trap regions within the Artemis exploration zone.We found that Amundsen, Cabeus, and Cabeus B craters retained the greatest quantities of ice potentially detectable with GPR.We found significant variance in model predictions for near-surface ice deposits, indicating that shorter term processes dominate ice retention in the upper 10 m.Although our model is inconclusive for surface level deposits, Faustini, de Gerlache, and Amundsen craters retained the greatest quantities in the upper 6 m and may be better targets for instruments with <10 m depth sensitivity than other cold traps in this study.Of the modeled cold traps, Shackleton was least likely to retain subsurface ice due to its young formation age, proximity to the pole, and lack of preserving ejecta layers deposited after its formation.Model variance due to the precise sequencing of cold trap formation, ejecta deposition, and ice delivery events will be constrained by dating samples returned by upcoming missions from the Artemis program.We showed that basin ice and ejecta delivery play crucial roles in retention of ice at the lunar south pole.Buried ice deposits beneath lunar polar cold traps have likely been exposed to reworking by ballistic sedimentation and are thinner and less extensive than previously reported.

Figure 1 .
Figure 1.Sequence of events resulting in layered stratigraphy in lunar south pole craters.Grayed out stratigraphy columns indicate that the column/crater does not yet exist.Notice how the presence and thickness of layers change with time as craters are formed and geologic events occur.We do not include ice modification processes in this figure or in the MoonPIES model.

Figure 2 .
Figure 2. Absolute model ages of craters and basins (inset) incorporated within the MoonPIES model.The ages of polar craters were drawn from Tye et al. (2015), Deutsch et al. (2020), and Cannon et al. (2020).Lunar epochs were defined by the ages of the basins(Orgel et al., 2018).Asterisks indicate a basin that deposits ejecta into at least one of the modeled craters.
and does not track ice once it leaves one of the model stratigraphy columns (see Sections 4 and 4.2).

Figure 4 .
Figure 4. Ballistic sedimentation volatilized fraction as a function of initial ejecta fraction (f ejecta ) and ejecta temperature.Volatilized fraction is expressed as the mean (left) and standard deviation (right) cold trap material exceeding 110 K in model simulations, computed over 50 runs with random initial ejecta positions.High temperatures and ejecta fractions result in high volatilized fractions, as expected.Ranges of volatilized fractions are indicated for ballistic sedimentation events resulting from craters (cyan) and basins (orange, present-cold-moon; red, warm-ancient-moon; Artemieva & Shuvalov, 2008) modeled in this work.

Figure 5 .
Figure 5. Cumulative Kinetic Energy (left) from each of the large age-dated craters, excluding basins.The kinetic energy contours are plotted over an average solar illumination map of the south pole (AVGVISIB_75S_120M_201608.LBL; Mazarico et al., 2011).Water ice cold trap extents (right; Landis et al., 2022).

Figure 7 .
Figure 7.A comparison of Faustini, Haworth and Cabeus cold traps over three different Monte Carlo model runs (including ballistic sedimentation effects).Shading indicates the quantity of ice in each layer (see reported values in FigureS7in Supporting Information S1).We note that the cold trap ages and stratigraphic sequence of a particular run differs due to random variation in ejecta ages and ice delivery.Runs A and B show typical columns for all three cold traps.In run C, the large ice layer in Faustini and Haworth is the result of a basin-scale cometary impact, while the absence of layering in Cabeus indicates that it formed later than all possible ejecta sources.

Figure 8 .
Figure 8. Model stratigraphy columns for the same model run with (a) ballistic sedimentation and (b) no ballistic sedimentation.Shading indicates the quantity of ice in each layer (see reported values in FigureS8in Supporting Information S1).When ballistic sedimentation was accounted for, large pure ice layers were lost from the base of the oldest cold traps(Faustini, Haworth, and Shoemaker).Shallower layers retain a similar ice % in both cases.Although this is a representative outcome, it should be noted that the absolute quantity of ice and stratigraphic sequence changes from run to run (Figures7-9).

Figure 9 .
Figure9.Sum total equivalent ice thickness retained in each cold trap stratigraphic column across 10,000 runs grouped by formation era and sorted by latitude.The width of each violin is scaled by the total number of runs retaining at least 1 m of ice.Median and quartiles are indicated as dashed lines.Without ballistic sedimentation, total ice thickness is greater, particularly for Pre-Nectarian and Nectarian cold traps.The difference in ice retention is smaller for Imbrian cold traps where few basins and local impacts disturb ice.About 1,000 runs retained >1 m for Wiechert J and only 3 runs retained >1 m of ice in Shackleton.

Figure 10 .
Figure10.Depth and thickness of ice layers for 10,000 model runs with and without ballistic sedimentation (orange and blue, respectively), represented as kernel density estimation (KDE) contour plots.Histograms indicate the distribution of depths and layer ice thicknesses as counts without KDE smoothing applied.Shaded regions indicate "gigaton" zones, where layer thickness exceeds 10 m.For Cabeus, horizontal lines denote LCROSS maximum excavation depth(Luchsinger et al., 2021;Schultz et al., 2010).

Figure 11 .
Figure11.Boxplot of ice retained in the upper 6 m (top) and upper 100 m (bottom) of each modeled cold trap over 10,000 runs (for reference, the LCROSS impact excavated 6-10 m into Cabeus crater;Luchsinger et al., 2021;Schultz et al., 2010).The top of each box denotes the third quartile and whiskers denote the 99th percentile.Points denote outliers above the 99th percentile value.The dashed line indicates 5% ice concentration which was measured in the LCROSS ejecta plume(Colaprete  et al., 2010).