Polar Ice Accumulation from Volcanically Induced Transient Atmospheres on the Moon

Water ice exists at the lunar poles, but its origin, abundance, and distribution are not well understood. One potential source of water to the poles is the volcanic outgassing of volatiles from the lunar interior and subsequent condensation of erupted water vapor as surface ice. We investigate whether volcanic outgassing is a viable source for the accumulation of lunar polar water ice. We construct a model that accounts for volcanic outgassing, atmospheric escape to space, and surface ice accumulation over the period of peak lunar volcanic activity (4–2 Ga) and map the resulting water ice distribution and abundance using current surface temperature data from the Lunar Reconnaissance Orbiter. Our model suggests that ∼41% of the total H2O mass erupted over this period could have condensed as ice in the polar regions, with thicknesses up to several hundreds of meters. The south pole accumulates roughly twice the ice mass of the north, and the southern deposits are thicker. Typical modeled eruptions generate collisional atmospheres with lifetimes of ∼2500 yr. However, these atmospheres are episodic and generally do not persist between eruptions. Roughly 15% of an atmosphere’s water vapor mass forms a frost on the lunar nightside, while the transient atmosphere persists. Our work suggests that the volcanically active period of the early Moon would have been punctuated by short-lived, collisional atmospheres that enabled the efficient sequestration of large quantities (8.2 × 1015 kg) of water ice at the poles and the temporary diurnal availability of water ice and vapor at all latitudes.


Introduction
Water ice deposits exist at the lunar poles, where temperatures are low enough for these deposits to remain stable over long timescales. The Moon's current surface topography and low ∼1.5°obliquity produce permanently shadowed regions (PSRs) near the poles, the coldest of which (<110 K) are referred to as cold traps. The temperatures within these cold traps are low enough to allow water ice to remain stable for billions of years (Watson et al. 1961). Water ice has been both directly and indirectly observed in polar cold traps, but its overall abundance and distribution are uncertain (Colaprete et al. 2010;Mitrofanov et al. 2010;Fisher et al. 2017;Li et al. 2018;Rubanenko et al. 2019). Moreover, the origins of lunar ice are unclear and may have been a combination of several different sources. The most widely discussed mechanisms of volatile delivery are from asteroid and/or cometary impacts, solar wind ion implantation, and volcanic outgassing from the lunar interior (Arnold 1979;Morgan & Shemansky 1991;Prem et al. 2015;Lawrence 2017), each of which would affect the abundance, distribution, and composition of ice in unique ways (Lawrence 2017). Interest in the science and resource potential of lunar ice, especially in the context of NASA's upcoming Artemis missions, is driving efforts to further characterize volatiles and understand their history on the Moon.
We explore lunar volcanism as a source for the Moon's polar water ice deposits. Estimates of the volatile content of lunar magma suggest that water vapor was a significant component of the gases released during eruptions (Saal et al. 2008;Hauri et al. 2011;Robinson & Taylor 2014;Chen et al. 2015;Rutherford et al. 2017;Ni et al. 2019). In order to accumulate the water released in an eruption as ice, the total ice accumulation rate into polar cold traps must exceed the rate of water vapor loss to space. Recent work has discussed the possibility of ancient, volcanically induced, transient atmospheres, which may have played a roll in the delivery of water to the poles. Needham & Kring (2017) proposed that volcanism on the ancient Moon could have released enough volatiles to produce a fully collisional atmosphere at ∼3.5 Ga that would have dissipated over a period of ∼70 Myr. This regime would be entirely distinct from the surface-bounded water exospheres thought to interact with the Moon's cold traps in the present day . Aleinov et al. (2019) further characterized this atmosphere using a general circulation model to explore its viability and climate. Tucker et al. (2021) investigated the escape processes relevant to such an atmosphere and their effects on atmospheric lifetimes. Subsequent work by Head et al. (2020) suggested that the total volume of magma erupted was probably lower than that estimated by Needham & Kring (2017) and was released by smaller eruptions over a longer period of time from 4 to 2 Ga.
In this work, we model the accumulation of volcanically sourced ice in lunar polar cold traps, using an eruption timeline consistent with that proposed by Head et al. (2020) and an atmospheric escape model that considers Jeans, photodissociative, and sputtering escape. We compare the timescales required to condense water ice out of volcanic atmospheres with rates of water loss to space. We use bolometric temperatures derived from the Diviner instrument (Paige et al. 2010b) to determine where on the lunar surface water ice should condense and how much ice this mechanism could have delivered to the polar cold traps. This work evaluates the feasibility of volcanism as a source for lunar polar volatiles, characterizes the distribution of ice resulting from such a scenario, and demonstrates where thick ice deposits are most likely to be found on the Moon. Figure 1 shows the primary components of the model used in this work. The model balances volcanic outgassing with atmospheric escape to space and surface ice accumulation from 4 to 2 Ga. The model components are described in detail in the following subsections.

Volcanic Eruptions
The model simulates lunar volcanic eruptions during the Moon's most volcanically active period from 4 to 2 Ga based on the work by Head et al. (2020). The timing of each eruption is determined by a probability distribution derived from previous work on lunar volcanic activity (Head & Wilson 1992;Whitten & Head 2015), which is characterized by a peak occurring around 3.7-3.5 Ga followed by a rapid fall-off toward 2 Ga (Figure 2). The amount of gas released in an individual eruption is determined by the volume of magma erupted onto the surface and the volatile content of the magma.  predicted a range of erupted volumes of ∼10 1 -10 3 km 3 per eruption. The range of volumes estimated from observed lava flows agrees with this, with typical mean erupted volumes estimated to be ∼10 2 km 3 (Head 1976;Yingst & Head 1997, 1998Whitten et al. 2011). The volatile mass fraction of lunar magma has been estimated to range from 700 ppm, based on the size of lava ponds for sinuous rilleforming eruptions , to ∼3400 ppm, based on analysis of picritic magma sampled by Apollo missions (Rutherford et al. 2017). Here we use an average value of 2000 ppm for the volatile mass fraction after Head et al. (2020). This range of eruption volumes and volatile mass fraction leads to estimates of ∼10 10 -10 13 kg of gas released per eruption. The volatile content of each modeled eruption is set by the volatile fractions of lunar picritic magma: 41% CO, 33% H 2 O, 9.6% SO 2 , 4.9% H 2 S, 9.6% COS, and 1.5% F .
While there is considerable uncertainty in estimating the preeruptive water abundance within lunar magmas (Chen et al. 2015), this model's assumed H 2 O volatile fraction of 667 ppm falls within published estimates made using a variety of techniques (e.g., 745 ppm, Saal et al. 2008;615-1410615- ppm, Hauri et al. 2011110 ppm, Chen et al. 2015;801-923 ppm, Rutherford et al. 2017;84 ppm, Ni et al. 2019). The Needham & Kring (2017) work on volcanically induced atmospheres used a water fraction of 1.98-9.9 ppm, derived from the review of the topic made by Robinson & Taylor (2014). However, Robinson & Taylor (2014)   667 ppm H 2 O is therefore a reasonable assumption. Further discussion of model behavior under different assumptions of water abundances is included in Section 4. Head et al. (2020) used estimates of mare basin lava fill depths to calculate a total erupted lunar magma volume of ∼10 7 km 3 . Using a typical density for lunar basaltic magma of 3000 kg m −3 and a volatile mass fraction of 2000 ppm, the total mass of gas released by all lunar maria-forming eruptions is calculated to be ∼6 × 10 16 kg. Adopting this total mass, we explore model behavior under two different assumed sizefrequency distributions (SFDs) for eruption mass ( Figure 3): (1) an exponential distribution, P(m) ∝ (e − m/ β )/β, and (2) a  . Erupted vapor mass histograms for two different model runs. The blue circles are derived from an exponential SFD with a total of 50,000 eruptions. The orange triangles are derived from a power-law SFD with a total of 545,000 eruptions. The exponential distribution is the default assumption of our model. power-law distribution, P(m) ∝ m − α , where m is the erupted vapor mass, β is the mean erupted mass for the exponential distribution, and α is the power-law exponent. The SFD from which eruption sizes are pulled remains the same over a model run, though the number of eruptions in any time interval is determined by the time distribution shown in Figure 2. The exponential distribution uses a mean of β = 1.2 × 10 12 kg, which is consistent with typical erupted magma volumes. The long tail of the exponential distribution also allows larger, less frequent eruptions to occur . The constraint on the total vapor mass released over 4-2 Ga requires 50,000 eruptions to occur over this period, within the range of 30,000-100,000 eruptions estimated by Head et al. (2020). Alternatively, a power-law distribution (Figure 3) is similar to eruption magnitude-frequency relationships observed on Earth (e.g., Pyle 1995). In this case, we limit the smallest eruption volume to 10 km 3 , corresponding to a vapor mass of 6 × 10 10 kg, and tune the power-law exponent so that the distribution approaches zero near a mass of ∼10 13 kg, corresponding to the largest known lunar eruption volume (∼10 3 km 3 , Schroeter's Valley; Head et al. 2020). The resulting power-law exponent is α = 3.2, and the number of eruptions required to account for the total released vapor mass is 545,000. By exploring these two vapor mass probability distributions, we are able to probe model behavior when volcanism is characterized both by less frequent, larger eruptions and more frequent, smaller eruptions. Because the exponential distribution is consistent with the estimated total number of eruptions ) and the mean erupted magma volumes (Head 1976;Yingst & Head 1997, 1998Whitten et al. 2011;Head et al. 2020), we use the exponential distribution as the default mode of the model. The results presented hereafter were obtained using the exponential distribution, and the effects of using the power-law distribution are discussed in Section 4.

Atmospheric Escape
In competition with polar ice accumulation, the primary escape mechanisms for an ancient lunar atmosphere are thermal (Jeans) escape, photochemical loss, and solar wind sputtering. Each of these mechanisms is sensitive to atmospheric attributes such as total mass, temperature structure, composition, mean molecular mass, and chemistry, as well as external factors like the presence of a lunar magnetic field or the strength of the solar wind (Tucker et al. 2021). Some of the previous work concerned with atmospheric lifetimes of volcanically induced lunar atmospheres (Needham & Kring 2017;Head et al. 2020) used a constant escape rate of 10 kg s −1 , based on Jeans escape estimates made by Vondrak (1974). However, a constant escape rate does not reflect the sensitivity of escape mechanisms to small differences in atmospheric parameters that may differ from one eruption to another or change over time (Tucker et al. 2021). Therefore, we employ an escape model that allows for a range of escape rates that depend on atmospheric total mass, mean molecular mass, and temperature. We summarize each of these below.
The atmospheric escape model considers Jeans escape, photodissociative escape, and solar wind sputtering; details of this model are described in the appendices. The isothermal Jeans escape model assumes that the atmosphere is composed of a single species with a molecular mass that is the weighted average of all atmospheric species (Jeans 1921). Escape rates are determined at the exobase, the height of which depends on total atmospheric mass (see Appendix A). Exobase height affects both the surface area of the sphere over which escape can occur and the number density of particles available for escape. This height also affects the gravitational binding energies of the atmospheric particles subject to escape. Although the results presented below demonstrate that Jeans escape is not a dominant atmospheric escape mechanism (see Figure 11 in Appendix A), it is included in the model for completeness. The Jeans escape model does not consider individual escape rates for different atmospheric species but uses the mean molecular mass of the atmosphere to calculate a net escape rate, consistent with an assumption of a well-mixed atmosphere. The mean molecular mass of an atmosphere produced after a single eruption in our model is 31.4 g mol −1 (Rutherford et al. 2017;Head et al. 2020). However, the model tracks atmospheric total mass, H 2 O mass, and CO mass and adjusts the mean molecular mass as H 2 O or CO is lost to photodissociation or H 2 O ice formation. The model assumes an isothermal atmosphere with a temperature determined by radiative equilibrium with 25% lower solar irradiance at 3.5 Ga (Equation (A2)). This radiative equilibrium temperature is 253 K, which also falls between average temperatures determined by Aleinov et al. (2019) for wet and dry CO-dominated atmospheres around the Moon.
Our model of photodissociative escape follows the approach of Tucker et al. (2021). The model accounts for photodissociative escape of both CO and H 2 O. We estimate the photodissociation rate of H 2 O due to solar ultraviolet photons by examining a hypothetical lunar atmosphere composed entirely of H 2 O having a mass equal to that of the total water vapor in the volcanically induced atmosphere at each time step in the model (Figure 4). This approach assumes that the opacity of the other atmospheric species can be neglected at the relevant wavelengths. The photodissociation rate at a particular height is calculated as a function of the line-of-sight optical depth, H 2 O photodissociation cross section, atmospheric particle number density, and solar ultraviolet photon flux integrated over relevant wavelengths for H 2 O (Equations (A9) and (A11); Tucker et al. 2021). To calculate the total photodissociation rate on the lunar dayside, we use a constant exobase temperature and integrate the height-dependent volumetric rate upward from the exobase. Then, we multiply this rate by the hemispheric area at the exobase altitude to arrive at a total rate in molecules per second. We assume that 50% of photodissociated molecules escape, similar to Tucker et al. (2021). Note that this method does not account for the reduction of photodissociation rates due to the diffusion-limited transport of water vapor through the CO-dominated atmosphere, which would tend to raise the amount of ice accumulation at the poles. The model calculates photodissociation-driven mass-loss rates of CO in the same manner, though using photodissociation cross sections from Heays et al. (2017) and the TIMED SEE (Thermosphere Ionosphere Mesosphere Energetics and Dynamics Solar EUV Experiment) solar irradiance spectrum (Woods et al. 2005) for solar ultraviolet photon flux at photodissociative wavelengths for CO.
Atmospheric mass loss due to sputtering occurs when solar wind particles transfer momentum to atmospheric molecules, imparting them with escape energies and trajectories. Solar wind sputtering yields are defined as the number of atmospheric particles ejected per incident solar wind particle and are given for H + and He + solar wind ions in Table 4 of Tucker et al. (2021). The model assumes an early solar wind flux of 1.65 × 10 13 H + m −2 s −1 and that a fraction of 33% of the solar wind flux reaches the surface, consistent with a 0.5 μT average magnetic field strength for the early Moon (Garrick-Bethell et al. 2019). Sputtering rates depend on the crosssectional area of the lunar atmosphere in the solar wind as defined by the lunar exobase height and are therefore a function of total atmospheric mass (Equation (A13)). As atmospheric mass decreases, the exobase height decreases, causing sputtering rates to also decrease ( Figure 4). The CO photodissociative and sputtering escape rates calculated by Tucker et al. (2021) are 200 and 23 kg s −1 , respectively, for an atmospheric mass of 5 × 10 15 and temperature of 300 K. These rates are significantly higher than the range of rates calculated in this model, primarily due to the higher atmospheric mass and temperature used by Tucker et al. (2021).

Ice Accumulation
Surface temperature and atmospheric water vapor pressure are the two primary factors that control the condensation of ice onto the lunar surface. Ice will condense if the atmospheric water vapor density (ρ v ) is greater than the saturation vapor density (ρ sat ) at a surface, and vice versa for sublimation. ρ sat is a nonlinear function of temperature, so a small increase in temperature results in a large increase in ρ sat . Therefore, the maximum temperature of a region determines the stability of ice in that region (Watson et al. 1961;Paige et al. 2010b). We can calculate the rate of ice accumulation in a region with maximum temperature T max by considering a laminar air layer at the surface through which water vapor can diffuse and above which the atmosphere is well mixed (Wallace & Sagan 1979;Bapst et al. 2018;Wilcoski & Hayne 2020). The efficiency of water vapor transport between the atmosphere and the surface will be limited by the thickness d of the laminar layer and the molecular diffusivity D of H 2 O through the CO-dominated laminar layer. The mass flux of water vapor through the laminar layer, and therefore the ice accumulation rate, is expressed as where ∂m ice /∂t has units of kg m −2 s −1 and is positive when ice is accumulated. The laminar layer thickness d depends on the kinematic viscosity of CO. Both kinematic viscosity and diffusivity D are functions of temperature and pressure at the surface (Chen & Othmer 1962;Chittenden et al. 2008). The laminar layer thickness d also depends on an aerodynamic surface roughness parameter z 0 and wind speeds above the surface (Equations (B1) and (B4); Wallace & Sagan 1979;Bapst et al. 2018). The aerodynamic roughness parameter z 0 of a surface is usually estimated from measurements of the wind profile above the surface. In the absence of any such measurements for the early Moon, we use Mars-like wind speeds consistent with those measured by the Phoenix Mars Lander (4.7 m s −1 at 2 m height) and the surface roughness parameter z 0 = 5 mm calculated for the Phoenix landing site. Additionally, Hébrard et al. (2012) developed a method of calculating z 0 from rock abundance data using observations from both Mars and arid regions of Earth. Using the range of rock abundances measured at the Chang'E-3 landing site by the Yutu rover (Di et al. 2016) and the method of Hébrard et al. (2012), we estimate a range of z 0 values from 2.0 × 10 −3 m to 1.1 × 10 −1 m for the lunar surface at the Chang'E-3 site (see Appendix B for details). The median z 0 value for the Chang'E-3 site also agrees well with the z 0 calculated for the Phoenix landing site. We use this range of z 0 values as an estimate of the lunar variability of z 0 and investigate model behavior over this range. Note that our estimates of the laminar layer thickness d are likely overestimates due to the instabilities that would exist between . Atmospheric escape rate as a function of atmospheric mass for H 2 O (solid blue line) and CO (dashed-dotted orange line) photodissociative escape and sputtering escape (green dashed line). For photodissociation, the atmospheric mass represents the mass of that single species (H 2 O or CO) present in atmosphere. Jeans escape rates (not shown here; see Appendix A) vary from 10 −11 to 10 −5 kg s −1 over the range of atmospheric masses from 10 8 to 10 14 kg and are therefore negligible compared to photodissociative and sputtering escape rates. parcels of air above surface elements with drastically different temperatures in a thin lunar atmosphere.

Diviner Data and Ice Stability
As an input, the model uses maximum temperature maps from the Diviner Lunar Radiometer Experiment on board the Lunar Reconnaissance Orbiter (Paige et al. 2010a). The Diviner maps show bolometric temperatures calculated from mean radiances from 6 to 18 hr local time and 18-6 hr angle for polar summers, both in increments of degree of effective lunar solar longitude for 10 Draconic years. These maximum temperatures are in good agreement with those found in Williams et al. (2019). Each pixel has a spatial resolution of 0.01°in polar azimuthal equidistant projection, or ∼250 m. For further details on how these temperatures were derived and for the archived data, see Landis et al. (2022). Each pixel of the maps is treated as a potential ice reservoir with a cold-trapping efficiency determined by its maximum temperature. Diviner temperatures are adjusted to account for ∼25% lower solar luminosity at ∼3.5 Ga (Feulner 2012) by assuming that temperature scales with the fraction of present-day insolation to the 1/4 power ( Figure 5). This temperature adjustment increases the area available for cold trapping at the poles by about 27%. Inherent in the use of Diviner maps is the assumption that current lunar surface topography and obliquity are similar to those of the Moon's past. The implications of this assumption are discussed further in Section 4.
Ice accumulation at each potential surface reservoir (map pixel) is coupled to the modeled atmosphere using Equation (1). In each time step, the ice accumulation rate is calculated at each surface reservoir based on that reservoir's temperature and the local partial pressure of water vapor of the atmosphere, and then the ice mass in each reservoir is updated. As ice accumulates, the equivalent water vapor mass is removed from the atmosphere, and the mean molecular mass of the atmosphere is updated. The atmosphere allows communication between surface reservoirs after ice deposition, such that warmer cold traps can lose ice gradually to their colder counterparts. This communication occurs until the atmosphere ceases to be collisional, after which point ice is no longer added to or removed from surface reservoirs and any water vapor remaining in the atmosphere continues to escape to space. An atmosphere remains collisional as long as the mean free path of a particle at the surface is less than the density scale height of the atmosphere at the surface. Halting ice accumulation after an atmosphere ceases to be collisional has no significant impact on the total amount of ice accumulated in cold traps, because by then the atmospheric water content has been depleted.
We do not actively model ice sublimation to space after the atmosphere has escaped during the model run, which allows this model to provide estimates of the total ice supply to the poles without incorporating a more complex model for ice burial/destruction post-emplacement. However, after a model run has completed, we calculate the ice loss that would have occurred owing to sublimation if the ice were exposed to space for 4 Gyr (Hayne & Aharonson 2015) and estimate the amount of ice that would remain. The model time step is 2 yr, which is chosen to be roughly two orders of magnitude faster than water is removed from the atmosphere. Head et al. (2020) estimated that the duration of even the largest known lunar eruptions would have been less than a year. Therefore, we add all the vapor mass from an individual eruption to the atmosphere in a single time step and assume that the atmosphere is well mixed over this period. Unlike the work by Prem et al. (2015) and Aleinov et al. (2019), we do not explicitly model the dynamics of the atmosphere, but instead treat the atmosphere as one wellmixed reservoir that interacts with all surface ice reservoirs simultaneously. A Hadley cell with a high latitudinal extent is expected for a slow rotator like the Moon (Held & Hou 1980), and assuming Mars-like surface wind speeds (∼1 m s −1 ), Figure 5. Maximum surface temperatures measured by the Diviner Lunar Radiometer Experiment and adjusted to account for ∼25% lower solar luminosity at ∼3.5 Ga from ±60°latitude to the poles in the south (left) and the north (right). efficient atmospheric mixing would occur on timescales much less than a year.
An additional effect we consider is the reduction of the total atmospheric water vapor content due to the condensation of water frost onto the lunar nightside. Before running the full 2 Gyr model, we run a series of high temporal resolution ice accumulation models to determine what fraction of the atmospheric water vapor condenses on the nightside given an erupted water vapor mass. We run these models with a time step much shorter than one lunar day and allow the temperature at each point on the lunar surface to change with time using the analytic temperature function derived by Hurley et al. (2015) adjusted for 25% lower solar luminosity. We assume that the presence of an atmosphere does not have an effect on surface temperatures, since CO is radiatively neutral (Aleinov et al. 2019), and we assume that the sensible heat of the atmosphere is negligible owing to its low density. We let each of these models equilibrate for several lunar days and then record the fraction of the total atmospheric water vapor mass condensed as ice on the nightside. These models are run for a range of total atmospheric masses from 10 6 to 10 16 kg. We use the results of these high temporal resolution models to determine the fraction of water vapor sequestered on the nightside at a particular time for any total atmospheric water vapor mass. In the full 2 Gyr model run, the available atmospheric water vapor mass that is able to interact with cold traps (i.e., not trapped as frost on the nightside) is adjusted at each time step based on the total atmospheric water vapor mass. For instance, if an atmosphere has a total water vapor mass of 10 11 kg but 15% of this mass is condensed as frost on the nightside, then the ice accumulation rates in polar cold traps are calculated using an available atmospheric water vapor mass of 0.85 × 10 11 kg. The results and implications of the nightside condensation model are discussed further in Sections 3 and 4.

Atmospheric Lifetime and Ice Condensation Timescale
Figure 6(a) shows the total mass and H 2 O mass of the lunar atmosphere over 4-2 Ga for a typical model run. Each spike in total atmospheric mass is an eruption occurring at that time step. Concurrent spikes of H 2 O occur based on its erupted mass fraction. Typical atmospheric masses produced during an eruption are ∼1.2 × 10 12 kg (∼0.5 μbar surface pressure), corresponding to the mean erupted vapor mass consistent with Head et al. (2020). This typical atmospheric mass is four orders of magnitude lower than that of the atmosphere proposed by Needham & Kring (2017) and subsequently explored by Rate of water removal from the atmosphere due to ice accumulation (purple) and photodissociation (yellow) over time. As shown in the insets, when atmospheric water vapor mass is high, ice accumulation dominates over photodissociation and vice versa. Aleinov et al. (2019) and Tucker et al. (2021). After an eruption, the atmosphere begins to lose mass owing to escape to space and accumulation of water ice at the poles. Total atmospheric escape rates are the sum of escape rates due to Jeans, photodissociative, and sputtering escape. Each of these processes is a function of total atmospheric mass, which also controls the exobase height. Atmospheres remain collisional at the surface until they fall below ∼10 7 kg in total mass. Atmospheric escape is dominated by photodissociation and sputtering ( Figure 4).
The lifetime of an atmosphere with the mean erupted mass of 1.2 × 10 12 kg is ∼2500 yr. This is shorter than typical intervals between eruptions (∼22,000 yr), even at ∼3.5 Ga, when eruptions are more frequent. The atmospheres produced by modeled eruptions are episodic and do not represent a continuous collisional atmosphere over longer periods, though numerous collisional atmospheres appear for ∼1000 yr periods. The timescale for H 2 O to be removed from a typical atmosphere due to photodissociative escape to space and condensation of ice at the poles is about 50 yr.
The total vapor mass released by eruptions from 4 to 2 Ga is 6.0 × 10 16 kg ; of this, 2.0 × 10 16 kg is water vapor based on our assumptions of erupted volatile mass fractions . H 2 O photodissociation rates are comparable to net ice accumulation rates, so the total mass of ice accumulated will depend on the balance between these two processes. Figure 6(b) shows the rate of water removal from the atmosphere due to both photodissociation and ice accumulation. When erupted H 2 O masses are larger, ice accumulation dominates over photodissociation, and vice versa. Therefore, when larger eruptions dominate (exponential SFD), a significant fraction of ice accumulates (∼41% of total H 2 O mass). Because ice accumulation rates depend on the assumed aerodynamic surface roughness z 0 , using our calculated range of z 0 values gives a range of accumulated ice fractions of 39%-52%. H 2 O photodissociation rates may be overestimated here because we do not consider the finite timescale of transport of H 2 O from the lower to the upper atmosphere, which could in reality limit the supply of H 2 O available for photodissociative escape. Along with our use of maximum temperatures to calculate ice accumulation rates, this means that these values can be considered lower bounds on the total fraction of erupted H 2 O that can be accumulated as ice in polar cold traps using an erupted H 2 O volatile fraction of 667 ppm for lunar magma. The total amount of ice mass erupted over the 4-2 Ga period depends on the assumed water fraction outgassed by lunar magma. As the outgassed water fraction decreases, photodissociation rates will dominate over ice accumulation rates, and a smaller fraction of erupted water vapor will accumulate as ice.
Results of the nightside condensation model indicate that roughly 15% of the total atmospheric water vapor mass should be condensed as ice on the lunar nightside at any given time. This condensed fraction varies slightly from ∼12% to ∼17% for total atmospheric water vapor masses from ∼10 6 to ∼10 12 kg, respectively. The fraction of atmospheric water vapor that is able to condense on the nightside is limited by the rate at which water vapor can diffuse through the laminar layer to the surface. Figure 7 shows the thickness of this frost as it varies with latitude and time for an arbitrary longitude and an initial atmospheric water vapor mass of 4.0 × 10 11 kg. A frost is formed at all latitudes on the nightside owing to the Moon's large diurnal temperature swings, with temperatures reaching as low as 95 K at the equator (Vasavada et al. 2012). Frost begins to accumulate after sunset and increases in thickness throughout the night, becoming concentrated at the dawn  terminator. The thickness of frost at the equator reaches 5 μm assuming an ice density of 920 kg m −3 , before beginning to ablate at dawn. The latent heat absorbed by the surface owing to ice accumulation would lead to a negligible surface temperature change of <0.1 K and is therefore not accounted for in the nightside condensation model. Note that while the frost is thicker at higher latitudes, the majority of the water mass exchange occurs at lower latitudes owing to their greater surface area. At latitudes poleward of about ±80°the frost forms caps that do not ablate away during the day. Eventually, these caps will recede as ice accumulates in permanent cold traps and the water vapor mass of the atmosphere decreases, causing ice to become unstable at these latitudes on surfaces illuminated during the day. Equatorial frost reaches roughly 16°in longitude (∼500 km) past the terminator on the lunar dayside before being completely ablated. The frost would be optically thick, and would therefore create a bright deposit extending outward from the dawn terminator, along with visible frost caps at high latitudes. The reduced mass of atmospheric water vapor due to nightside condensation decreases ice accumulation rates in permanent cold traps but does not ultimately prevent a significant amount of water from reaching polar cold traps. Figure 8 shows the distribution and thickness of ice deposits resulting from volcanic outgassing over 4-2 Ga assuming the exponential eruption SFD. The thickest deposits are on the order of hundreds of meters thick, assuming an ice density of 920 kg m −3 . Most of the ice lies poleward of 80°latitude at both poles (Figures 8(c) and 8(d)). Ice thickness is not uniform among ice reservoirs, and ice is concentrated in the coldest cold traps, which accumulate ice the fastest. This occurs because the accumulation rate depends on the diffusivity, laminar layer thickness, and saturation vapor pressure, all of which are complex functions of temperature. Additionally, warmer ice reservoirs tend to gradually lose ice to colder ice reservoirs via the atmosphere, so the longer a water-bearing atmosphere persists, the more concentrated ice becomes. The model does not consider processes that destroy or bury ice (e.g., impact gardening; Arnold 1979) once the atmosphere becomes noncollisional. Therefore, ice thicknesses represent the amounts of ice that could have been supplied by volcanic outgassing, and not necessarily the net ice amounts that still exist in the present day.

Distribution and Abundance of Polar Ice
More ice accumulates at the south pole than the north. This is due to both the greater total cold trap area in the south and the greater area of cold traps with very low temperatures (Williams et al. 2019;Hayne et al. 2021). The south polar region accumulates ∼1.7 times the ice mass of the north, with 5.19 × 10 15 kg (83% poleward of 80°) and 3.04 × 10 15 kg (86% poleward of 80°) in the south and north, respectively. Ice deposits cover areas of 1.02 × 10 5 km 2 (72% poleward of 80°) and 8.16 × 10 4 km 2 (78% poleward of 80°) in the south and north, respectively. The mean and maximum thicknesses of ice deposits are 49 and 410 m, respectively. The south polar ice deposits also tend to be thicker than the north polar deposits because there are more cold traps with extremely low temperatures in the south. These deposits are thickest primarily because these cold traps accumulate ice the fastest, and also because they area able to scavenge ice from warmer cold traps while a collisional atmosphere allows communication of vapor between cold traps. The thickest ice deposit predicted by the model lies within an exceptionally cold crater that itself lies within Haworth crater (87.5°S, 347.6°E). The mean and maximum thicknesses of south polar ice deposits are 55 m (63 m poleward of 80°) and 410 m, respectively. The mean and maximum thicknesses of north polar ice deposits are 41 m (45 m poleward of 80°) and 319 m, respectively. Figures 8(e) and 8(f) show ice deposits remaining after 4 Gyr of sublimation to space, assuming that ice is continuously exposed at the surface, with no burial by crater ejecta and no formation of sublimation lag; this is intended to capture thermal effects on exposed ice only. After long-term sublimation, ∼2.2 times more ice remains in the south (3.57 × 10 15 kg) than the north (1.65 × 10 15 kg), and the areas covered by ice are reduced by factors of 4.1 and 6.7 in the south and north, respectively. The mean thickness of the remaining ice deposits is 152 m, while the maximum thickness is unaffected by sublimation to space. Figure 9 shows the fractional surface area covered by ice deposits as a function of latitude at both poles. The fractional surface area of ice follows a strong increasing trend with proximity to the pole in the south. However, the poleward trend is much more subdued in the north, with a peak of ice concentration shifted to ∼85°latitude. This asymmetry between the two poles appears even when sublimation to space over 4 Gyr is taken into account.

Discussion and Conclusions
This model suggests that volcanism could have been a viable source for lunar polar water ice deposits. Individual eruptions very likely released enough vapor mass to have created collisional atmospheres around the Moon. The fraction of erupted H 2 O accumulated as polar ice in an individual eruption depends on the erupted H 2 O mass. For the exponential erupted vapor mass distribution consistent with Head et al. (2020), ∼41% of the total H 2 O mass erupted over the 4-2 Ga period is deposited as ice at the poles. Photodissociation rates in our model are likely overestimates, and therefore this accumulated ice fraction can be considered an underestimate. Photodissociative and sputtering escape processes dominate over Jeans escape, and atmospheric lifetimes resulting from typical eruptions are ∼2500 yr. Typical eruption intervals at peak volcanic activity (4-3 Ga) are ∼22,000 yr, so it is unlikely that a collisional atmosphere would have existed continuously over longer periods of the Moon's past. In any individual eruption, H 2 O is removed from the atmosphere owing to ice accumulation and photodissociation in about 50 yr. Even though atmospheres do not typically persist between eruptions, water is efficiently sequestered at the poles well before individual atmospheres escape.
Ice distribution and thickness are directly related to surface temperature because temperatures control ice accumulation rates. Therefore, ice is concentrated in the coldest cold traps (Figure 8). Ice thicknesses output by the model represent the potential supply of ice to cold traps due to volcanic outgassing. Ice thicknesses average 49 m after a typical model run, and 152 m if sublimation to space occurs over 4 Gyr. The south pole accumulates ∼2 times the ice mass of the north and retains a greater proportion of the remaining ice if long-term sublimation to space occurs following accumulation. Ice deposits at the south pole are thicker on average than those at the north, because the south polar cold traps tend to be colder. Therefore, the coldest cold traps in the south polar regions provide the best locations to search for volcanically sourced ice deposits.
The particular SFD of volcanic eruptions could have a meaningful impact on the accumulated polar ice. A power-law probability distribution for erupted vapor mass enables us to probe model behavior when smaller, more frequent eruptions dominate the total erupted vapor mass over the Moon's volcanically active period. When smaller eruptions dominate, ∼14% of the total erupted H 2 O mass is deposited as ice, compared to ∼41% for the exponential SFD. This smaller accumulated ice fraction occurs because smaller eruptions create atmospheres with lower water vapor masses. As atmospheric water vapor mass decreases, ice accumulation rates decrease more rapidly than H 2 O photodissociation rates, allowing photodissociative escape to dominate over ice accumulation for lower-mass atmospheres. The fraction of the total erupted H 2 O accumulated as ice in the power-law SFD ranges from 12% to 23% depending on the assumed aerodynamic surface roughness parameter z 0 . In the powerlaw case, the lower total mass of accumulated ice results in ice deposits with mean and maximum thicknesses of 11 and 146 m, respectively. The ratio of ice mass in the south to the north remains roughly the same in both eruption SFD models. The areal extent of ice deposits increases for the power-law SFD to 1.53 × 10 5 km 2 and 1.29 × 10 5 km 2 in the south and north, respectively, because smaller, shorter-lived atmospheres provide less time for ice to concentrate in the coldest cold traps. The total amount of volcanically sourced water ice that exists on the Moon today may therefore be sensitive to the dominant eruption sizes during peak lunar volcanism. Additionally, if the outgassed water fraction is lower (or higher), ice accumulation will behave more like it does in the power-law (or exponential) model case.
An important assumption used by the model is the estimated pre-eruptive water abundance of lunar melt. The water abundance controls the total amount of water released over the Moon's volcanic history, and therefore the amount of water available to be trapped at the poles. In addition, when a higher water content is released during an eruption, a larger fraction of that water is able to accumulate at the poles. While our assumption of 667 ppm by mass of water in lunar magma is reasonable, estimates of this value range from tens to thousands of ppm (Robinson & Taylor 2014). Previous work on ancient lunar atmospheres assumed a very low water content (Needham & Kring 2017;Aleinov et al. 2019;Tucker et al. 2021), and Renggli et al. (2017) note that hydrogen species other than water may need to be taken into account to avoid overestimating water content in lunar melts. We investigated model behavior under a range of assumed water abundance from 10 to 3500 ppm. Over this range of water abundances, the accumulated ice fraction, the total accumulated ice mass, and the mean ice thickness ranged from 5.8% to 59%, from 1.8 × 10 13 kg to 6.1 × 10 16 kg, and from 0.10 to 360 m, respectively, assuming the same abundances of nonwater volatiles. The higher of these extremes implies ice deposits significantly thicker than deposit thicknesses suggested by current observations (e.g., Rubanenko et al. 2019), while the lower extreme implies not enough volcanically sourced ice to survive ice-loss processes. If future observations determine the abundance of volcanically sourced ice within lunar volatile deposits, it could provide a means of constraining the erupted water content of lunar volcanic gases.
The fact that massive, exposed ice deposits are not observed on the Moon today means that if they exist, then they must be buried or mixed with regolith (e.g., Costello et al. 2021). Some layers within the modeled deposits could have been thick enough and accumulated rapidly enough to have remained relatively pure before another eruption occurred. The thickest (410 m) deposits would have accumulated 8.2 mm of ice per eruption, or once every ∼22,000 yr given the average eruption interval during the first 1 Gyr. Speyerer et al. (2016) predicted that the upper 2 cm of lunar regolith should be reworked in 81,000 yr, which is consistent with Costello et al. (2020). Over 81,000 yr, the coldest ice reservoirs should accumulate about 3.0 cm of ice, slightly greater than the depth predicted to be overturned by impact gardening over the same time period. This suggests that the stratigraphy of these deposits could generally be characterized by mixtures of ice with regolith, perhaps punctuated by pure ice layers resulting from large eruptions or clusters of eruptions. However, if any of these deposits have survived to the present day, their upper 3-10 m should have been intimately mixed with regolith owing to prolonged exposure to impact gardening (Costello et al. 2021).
While the only ice loss process modeled here is sublimation to space, other loss processes could have diminished the volcanic ice abundance retained at the poles. Calculations by Farrell et al. (2019) suggest that micrometeoroid impact ejection and vaporization of ice could be effective mechanisms of removing ice from polar cold traps. Assuming present-day impact flux and a surface composed of 100% ice, the rate of ice loss would be 3.38 × 10 −17 m s −1 based on the Farrell et al. (2019) estimates of H 2 O-loss rate. However, the impact flux at the Moon was likely ∼50-100 times higher at 3.9 Ga than it has been over the past 3.0 Gyr (Fassett & Minton 2013). Assuming that the impact flux was 100 times greater from 4.0 to 3.0 Ga, about 2.3 mm of ice would be removed every 22,000 yr, given a surface composed of 100% ice. This would effectively remove ice from average thickness deposits, which would accumulate ∼1.0 mm of ice on average every 22,000 yr. Deposits with a final modeled thickness greater than 115 m (e.g., Figures 8(e) and (f)) would still retain their ice. However, lunar impact flux dropped rapidly after ∼3.9 Ga, while volcanic activity on the Moon peaked at ∼3.6-3.5 Ga (Fassett & Minton 2013). Therefore, it is likely that ice loss due to micrometeoroid bombardment was most significant during the first several hundred million years of the modeled time period, before the majority of eruptions had occurred. Additionally, while these calculations assume a surface of pure ice, increased impact flux should also have increased rates of impact gardening and burial of ice from nearby ejecta, decreasing the amount of ice exposed at the surface and vulnerable to destruction by micrometeoroid bombardment. Therefore, it is possible that thicker ice deposits would have been able to survive loss due to these processes.
Although our model assumes that present-day surface temperatures are similar to those at the time of peak lunar volcanic activity, it is possible that the lunar obliquity was in fact higher during this period than the present-day ∼1.5°. It is also possible that changes in obliquity or true polar wander (TPW) could have altered the Moon's surface temperatures since volcanic volatiles were deposited (e.g., Siegler et al. 2011Siegler et al. , 2016. However, the timing and amplitudes of these dynamical changes are not well understood, and it is nonetheless likely that the Moon existed in a low-obliquity state (<10°) during most of the first billion years of its history (Siegler et al. 2015). Furthermore, much of the highlands terrain encompassing the Moon's polar regions was formed by ∼3.9 Ga, including large cold-trapping craters such as Haworth (∼4.18 Ga), Shoemaker (∼4.15 Ga), Faustini (∼4.10 Ga), and Amundsen (∼3.9 Ga) . Therefore, it is reasonable to assume that the present-day cold traps existed in more or less their present-day configuration prior to the Cassinistate transition occurring at roughly 2 Ga (Siegler et al. 2015). Because of the high rate of regolith overturn, any accumulated ice should have been buried or mixed with regolith and therefore protected from rising temperatures during subsequent obliquity excursions.
TPW could have slightly altered the locations of polar cold traps; however, the mechanisms of ice accumulation and atmospheric loss modeled here would have operated in the same way. There is considerable overlap between present-day cold traps and those predicted by Siegler et al. (2016) for a paleo-pole, meaning that much of the ice would accumulate in similar regions to those predicted by this model (see Siegler et al. 2016, Extended Data Figure 3). If volcanic outgassing occurred coincidentally with a TPW scenario, this model provides a mechanism that could help explain the observed neutron data as described by Siegler et al. (2016). The addition to this model of the temperature maps from a TPW scenario would provide estimates of the ice locations and thicknesses expected if such an event occurred. If future observations conclusively identify buried ice deposits, their presence could help constrain ancient lunar surface temperatures and possibly the location of a hypothetical paleo-pole.
Our results show that the fractional surface area of ice increases with a clear trend toward the south pole, but not so for the north. Intriguingly, Rubanenko et al. (2019) found a similar trend in crater depth-to-diameter ratio (d/D), with mean d/D decreasing toward the lunar south pole but not the north. In addition, they found that pole-facing slopes of craters tended to be slightly shallower than equator-facing slopes in the south. They interpreted these observations as evidence for infill of south polar craters with buried ice deposits. The results of our model are consistent with this interpretation. Rubanenko et al. (2019) found south polar craters with d/D on average ∼10 m shallower than those at lower latitudes. If the shallowing of these craters is due to ancient ice sourced from lunar volcanism, the ∼50 m ice thicknesses predicted by our model would be more than enough to account for these observations. This would imply that subsequent alteration, removal, and/or destruction of ice ablated volcanically sourced ice deposits to ∼20% of their original thicknesses. Alternatively, our assumed water fraction of erupted volatiles could be too high, or higher obliquity could have led to higher shadow temperatures, leading to overestimates of the total amount of water released or accumulated. Finally, this model predicts the accumulation of ice over significant areas where it would only be temporarily stable. Each eruption is followed by periods where water frost forms at all latitudes on the lunar nightside and persists even during daytime at high latitudes. Evidence for the past and present existence of H 2 O outside of cold traps, such as observations of hematite (Li et al. 2020) or the detection of molecular water (Honniball et al. 2021), suggest potentially widespread interactions between ice or water vapor and the lunar surface. If mechanisms of surface-water interaction on the sunlit Moon were operating in the Moon's past, similar to the mechanisms responsible for these modern observations, the increased water abundance due to frost may have increased the amount of surface-water interaction. These surface-water interactions could have created detectable chemical signatures (e.g., Li et al. 2020;Honniball et al. 2021) outside of cold traps or permanently shadowed regions.
Part of this work was supported by the NASA Lunar Reconnaissance Orbiter project.

Appendix A Atmospheric Escape Model
Each escape mechanism depends on the exobase height, which varies with total atmospheric mass. The exobase height is defined as the altitude above the surface at which the density scale height of the atmosphere is equal to the mean free path of an atmospheric particle (Pierrehumbert 2010). In other words, the exobase is the altitude above which the atmosphere is no longer collisional. The atmospheric scale height (H) varies with altitude above the surface and is convenient to describe as a function of radius from the Moon's center (r): where r 0 is the radius of the Moon, H 0 is the atmospheric scale height at the surface, R is the universal gas constant, T 0 is the temperature of the atmosphere at the surface, g 0 is the acceleration due to gravity at the surface, and μ is the molar mass of an atmospheric particle in units of kg mol −1 . The model assumes an isothermal atmosphere with a temperature (T 0 ) determined by radiative equilibrium with 25% lower solar irradiance at 3.5 Ga (Feulner 2012). The equilibrium temperature is given by 1 is the surface albedo (Vasavada et al. 2012), F sol is the solar constant, f = 0.75 is the fraction of present-day solar irradiance at 3.5 Ga, ò = 0.95 is the surface emissivity, and σ B is the Stefan-Boltzmann constant. The mean free path (λ mfp ) of the atmosphere depends on the number density of particles (n), which itself is a function of altitude above the surface. λ mfp and n are defined in Equations (A3) and (A4) as functions of radius from the center of the Moon (r): where σ coll is the effective collisional cross section of an atmospheric particle and n 0 is the particle number density at the surface. n 0 is determined by the total mass of the atmosphere (M atm ): where N A is the Avogadro constant. By equating the atmospheric scale height (Equation (A1)) and the mean free path (Equation (A3)), we can solve for the exobase height (or radius from the Moon's center) r ex : The result is a transcendental equation that we solve using Newton's root-finding method for a range of different atmospheric masses. Figure 10 shows exobase height as a function of atmospheric mass for isothermal atmospheres with T 0 = 253 K. Note that the exobase height for H 2 O is much higher than that for CO because λ mfp depends exponentially on the molar mass μ.

A.1. Jeans Escape
At each time step of a model run, the Jeans escape rate is calculated as a function of total atmospheric mass (M atm ), which controls the exobase height (r ex ) and therefore the number density of particles at the exobase (n(r ex ) = n ex ). Jeans escape occurs for atmospheric particles that have velocities higher than the gravitational escape velocity (Jeans 1921). Escape trajectories occur when these particles do not collide with any other atmospheric particles before escaping. Therefore, Jeans escape rates are calculated at the exobase. is the most probable velocity of a particle, k B is the Boltzmann constant, T ex is the temperature at the exobase, and λ esc is the Jeans escape parameter. λ esc is the ratio between gravitational binding energy and thermal energy for a particle at the exobase and is given by where G is the gravitational constant and M moon is the mass of the Moon. The final Jeans escape rate (kg s −1 ) is calculated by multiplying the Jeans escape flux (kg m −2 s −1 ) over the spherical surface area of the lunar exobase ( pr 4 ex 2 ). Figure 11 shows the resulting Jeans escape rate as a function of total atmospheric mass.

A.2. Photodissociative Escape
The photodissociative escape model is largely based on that of Tucker et al. (2021). The model estimates photodissociative escape rates for CO and H 2 O by considering hypothetical atmospheres composed solely of CO and H 2 O, respectively. The primary difference between these two hypothetical atmospheres will be the exobase height. The photodissociation rate at a given height above the surface depends on the line-ofsight optical depth (Tucker et al. 2021 where σ PA (λ) is the photoabsorption cross section of H 2 O or CO and is a function of photon wavelength λ, θ Z is the solar zenith angle taken to be 60°, and Chp is the Chapman grazing incidence integral (Smith & Smith 1972). The Chapman integral allows the line-of-sight optical depth to be calculated from the vertical column density N(r), which itself can be  The thickness d of the laminar layer can be estimated as where ν is the kinematic viscosity of the atmosphere, which depends on temperature and pressure above the surface (Chittenden et al. 2008). The molecular diffusivity D of water vapor in the laminar layer is calculated using the analytical solution for diffusivity in a binary gas system devised by Chen & Othmer (1962): The eddy diffusivity of water vapor in the turbulent layer is an order of magnitude larger than the molecular diffusivity in the laminar layer. Therefore, the ice accumulation rate is limited by the molecular diffusivity, and we can assume that the atmosphere is well mixed above the laminar layer.