Gradual Sequestration of Water at Lunar Polar Conditions due to Temperature Cycles

Migration of water molecules into the porous lunar soil can lead to sequestration of adsorbed water or ice at depth. Here, this process is modeled at several potential landing sites in the south polar region of the Moon. Desorption times are parameterized based on the Brunauer–Emmett–Teller (BET) isotherm model. Water sequestration is significant at depths where subsurface temperatures remain below the cold-trapping threshold, a condition satisfied at one of three sites considered, west of Haworth Crater. Model calculations predict a hydrated layer below a desiccated layer. The thickness of the desiccated layer is on the scale of the thermal skin depth, which for dust is typically centimeters. The underlying layer, decimeters thick, can be hydrated over a timescale of a billion years, reaching abundances on the order of 1 wt%. This sequestration process potentially simultaneously explains excess hydrogen concentrations outside of cold traps and the observed presence of a desiccated layer above a hydrogenous layer.


Introduction
The distribution and abundance of water in the polar regions of the Moon is of great interest for science and exploration but remains poorly understood (Lawrence 2017). Over time, water is delivered to the Moon by impactors, generated on the lunar surface by solar-wind interactions, or supplied by outgassing from the interior (Arnold 1979). Most of this water is quickly lost, but a fraction can be retained in the form of ice. In permanently shadowed cold traps, the sublimation rate of water ice into vacuum is negligible (Arnold 1979). Water might also be sequestered outside of cold traps through a process called "vapor pumping" or by sudden burial. Figure 1 provides a schematic overview of H 2 O sequestration mechanisms that have been proposed for the lunar polar regions. The topic of the present study is vapor pumping by temperature cycles (Figure 1(b)). If the subsurface is sufficiently cold, water molecules that randomly migrate into the porous subsurface are stored at depth. The rate of this sequestration is enhanced by diurnal temperature cycles; the periodically high temperatures allow the molecules to migrate downward. This sequestration mechanism works best with a quasi-continuous supply of water molecules to the surface, from delivery by a water exosphere or generated locally by the solar wind . Even low surface concentration amplitudes (less than a monolayer) suffice for vapor pumping, and some degree of diurnal variation of adsorbed water molecules, implying their emission, is suggested by observations (e.g., Hendrix et al. 2019). Thermal pumping can result in H 2 O concentrations far in excess of the surface concentration.
For airless bodies, this phenomenon was pointed out by Schorghofer & Taylor (2007), but it has been recognized earlier for Mars (Mellon & Jakosky 1993), where the source of the water is a humid atmosphere instead of adsorbed molecules on the surface. Schorghofer & Aharonson (2014) and Schorghofer & Williams (2020) published maps of where on the Moon this sequestration might occur. Recently, Reiss et al. (2021) also modeled vapor diffusion and adsorption, but not desublimation. None of the previous studies has obtained a vertical concentration profile for ice sequestered through this process from a steady supply of water molecules to the surface. Schorghofer & Taylor (2007) considered either constant temperature or a temporary ice cover; Schorghofer & Aharonson (2014) determined the yield of the process using a boundary-value formation. Siegler et al. (2015) considered pumping from an ice cover, and Reiss et al. (2021) obtained concentration profiles for adsorbed water but not for ice. Here, an improved model of subsurface vapor diffusion and adsorption is developed and applied to prospective landing sites to predict vertical concentration profiles and estimate expected absolute abundances of subsurface H 2 O outside of cold traps.

Improved Quantification of Desorption Rates
A quantification of molecular residence times or desorption rates that spans concentration regimes from a fraction of a monolayer to ice is necessary. The desorption or sublimation rate into vacuum is the number of molecules that leave a surface per area and time S = −d θ/dt, where θ is the areal adsorbate concentration. It can be written in terms of the equilibrium vapor pressure p as where k is the Boltzmann constant, m the mass of a H 2 O molecule, and T temperature. For ice, this is known as the Hertz-Knudsen equation (Persad & Ward 2016), but it applies more generally, as long as the velocity distribution of desorbed molecules from a solid wall is the same as for ice (Armand 1977;Oura et al. 2003).
A monolayer θ m ≈ 10 19 molecules m −2 . The specific surface area measured for the submillimeter soil fraction of Apollo samples is typically 500 m 2 kg −1 (Carrier et al. 1991). A monolayer of H 2 O weighs about 0.3 mg m −2 and hence represents a concentration of 150 ppm in weight. Cadenhead & Stetter (1974) measured a reversible isotherm for a lunar sample, p/p ice = f (θ/θ m ), where p ice is the equilibrium vapor pressure of ice. Schorghofer & Aharonson (2014) fitted this isotherm with a mathematical form that was empirical. Here, a physics-based isotherm model is used instead to fit this isotherm. Reiss et al. (2021) used the Langmuir adsorption model, where each adsorption site has the same binding energy, and once all sites are occupied no more molecules are adsorbed. This model is limited to a single monolayer and insufficient for the purpose of modeling ice sequestration.
For polar molecules and multilayer adsorption, the Brunauer-Emmett-Teller (BET) isotherm is appropriate: which has only one free parameter, c. For small p, θ/θ m → c p/p ice . For p close to p ice , θ/θ m → 1/(1 − p/p ice ). For one monolayer, ). Fitting of the measurements by Cadenhead & Stetter (1974) yields c = 8.231 (Figure 2). Because a measured isotherm is only available for a single temperature, the dependencies for T and θ are factored: is calculated with Equation (1) using vapor pressures from Murphy & Koop (2005).
The desorption rate S can be converted into a residence time τ. For less than a monolayer, τ = θ/S, and for more than a monolayer, τ = θ m /S. The relation is ambiguous at θ ≈ θ m because of the partial overlap of layers, and here it is assumed that a second monolayer starts to form only after the first monolayer is complete. This makes it straightforward to distinguish between mobile (exposed) and immobile (covered) molecules in model calculations.
With the BET model, the residence time approaches the values τ(θ = 0) = c τ ice and τ(θ = ∞ ) = τ ice . Water desorption rates on lunar soils have been measured by temperatureprogrammed desorption experiments (Poston et al. 2015;Jones et al. 2020). These measurements reveal sites of high binding energies over a fraction of the grain surface. Such sites will soon be permanently occupied, and the next water molecule arriving at the site will experience lower binding energy. In the current lunar polar environment, permanently adsorbed water will reduce the range of available binding energies. The reversible adsorption isotherm characterizes the population of water molecules that participate in the migration process.

Surface and Subsurface Temperatures at Prospective
Landing Sites Three study sites are selected in the south polar region, all prospective landing sites ( Table 1). The exact coordinates will not be known until after the landings, so the selected coordinates only represent an area in the vicinity outside of cold traps. Moreover, temperatures vary over rover traverses.
Surface temperatures are obtained from Diviner radiometer data binned into 96 solar longitudes (local time) bins and 2 ecliptic longitude bins (winter and summer) using the data Figure 1. Schematic illustration of three types of water ice storage in lunar soil. (a) Cold traps lie within permanent shadow (gray) and are so cold that sublimation into vacuum is less than the amount that was delivered to the cold trap over time. (b) Adsorbed water molecules migrate into the porous subsurface until they reach depths beyond the influence of diurnal temperature variations, where they remain indefinitely. (c) Ice that was quickly buried retreats far slower than ice exposed to vacuum. Figure 2. Fit of the BET isotherm to adsorption isotherm measurements. The vertical axis uses a combination of variables that turns the BET isotherm into a straight line. p 0 is the saturation vapor pressure of water. The BET Equation (2) has only one free parameter. The measurements with the highest partial pressures were omitted from the fit, but nevertheless θ → ∞ as p → p 0 . products of Williams et al. (2019). Winter and summer temperatures are averaged to obtain a time series with 96 steps per lunation. Table 1 includes the mean, minimum, and maximum temperatures at each of the three study sites. At the study site west of Haworth Crater, the time average is 113 K. For comparison, the sublimation rate of crystalline ice at 109 K is 100 kg m −2 Gyr −1 and at 114 K it is 1000 kg m −2 Gyr −1 (Murphy & Koop 2005). Hence, this site has a mean temperature close to the cold trap threshold, whereas the other two sites are distinctly warmer.
Subsurface temperatures are obtained by solving the onedimensional heat equation, using the measured Diviner temperatures as an upper boundary condition and a geothermal heat flux of 0.018 W m −2 as a lower boundary condition. The depth-dependent thermal properties assumed are those of Hayne et al. (2017). For these thermal properties, the diurnal skin depth is about 5 cm. Because the thermal model is driven by measured surface temperatures, the choice of thermal properties affects only the vertical temperature profile but not the surface temperature amplitude. Changes in thermal conductivity caused by the ice content are neglected.

Equilibrium Concentrations
The concept of thermal pumping can be introduced with a one-dimensional model of diffusion between grain surfaces. The desorption rate S is a function of temperature T and areal adsorbate concentration θ. For ice, S is the sublimation rate into vacuum, which only depends on T. The flux between two opposing surfaces with adsorbed water is q q = - In a continuum formulation, this difference can be approximated with J = λ ∂S/∂z, where λ is the spacing between the two surfaces and z the vertical coordinate. Denote the time average of a variable X over a period P by Since λ is constant in time, the time-averaged flux is given by Equilibrium refers to a dynamical balance between the number of volatile water molecules on the surface and in the subsurface. In equilibrium, the net flux vanishes, and the time average of S is the same at all depths and determined by its value on the surface: The surface population grows and shrinks over each lunation, and, outside of cold traps, it is periodically depleted. To evaluate á ñ S surface ( ), Schorghofer & Aharonson (2014) used a time-dependent model of the surface population, but outside of cold traps, where no net accumulation occurs on the surface, the time average of the desorption rate is the time average of the infall rate. Hence, á ñ S surface ( ) has a set value, namely the continuous component of the infall rate, e.g., 1000 kg m −2 Gyr −1 . The specific value is uncertain (see, e.g., Arnold 1979). This approach neglects the detrimental effect that space weathering has on the surface population of H 2 O, but given that this supply rate is equivalent to only one-quarter of a monolayer per lunation, most of the destructive incoming particles will not interact with the volatile H 2 O molecules. Ice forms at depths z where á ñ á ñ  S S Tzt surface , ice ( ) ( ( )). The left-hand side is assumed to be the infall rate, and the righthand side is approximately S T max t ice ( ( )). To this approximation, subsurface cold trapping occurs at depths that never exceed a specific temperature, and, if delivery is continuous, that temperature threshold is the same as for surface cold traps.
In the subsurface, the average q á ñ S T, ( ) in Equation (7) requires knowledge of θ(z, t). At depths where the temperature amplitude is negligible, θ will be constant after it has reached its equilibrium value. To calculate the equilibrium value of θ even at depths where temperature varies, it is assumed that the local adsorption dynamics follows an isostere. An adsorption isostere is the (T, p) combination that leaves θ constant, just as an adsorption isotherm is the (p, θ) combination at constant T. The justification is that only a limited amount of vapor can diffuse within a single temperature cycle and hence θ must be close to constant over a cycle. Once a stationary state is reached, θ should be constant even at depths where the temperature amplitude does not vanish. The isostere assumption is not valid on the very surface, and there is a "breathing" skin depth, thinner than the thermal skin depth, where θ varies periodically.
Assuming the adsorption dynamics follows isosteres, θ is constant with time (but not with depth) and given by For the parameterization chosen, the T dependence and θ dependence are factored, Equation (3), and therefore To determine θ eqlbr , á ñ S T ice ( ) is evaluated at each depth and then θ eqlbr is calculated from Equations (2) and (9) for a given value of á ñ S surface ( ). Figure 3 shows the subsurface temperatures and the number of monolayers adsorbed in equilibrium. Within the thermal skin depth there is little adsorption because of devolatization by periodically high temperatures. At the study site west of Nobile (Figures 3(a) and (b)), where temperatures are above the coldtrapping temperature, the concentration reaches only a fraction of a monolayer. At the study site west of Haworth (Figures 3(c) and (d)), the subsurface concentrations are much higher. For 100 kg m −2 Gyr −1 , the maximum concentration is about half a monolayer; for 1000 kg m −2 Gyr −1 , even ice is stable in the subsurface. (In equilibrium, all of the pore spaces would be filled with ice, but, as will be demonstrated below, there is not nearly enough time to reach this equilibrium.) At increasing depth, the geothermal gradient leads to decreasing adsorbate concentrations. The layer with significant hydration is decimeters thick. It starts below the influence of the diurnal thermal cycle and tapers out at greater depths due to the geothermal temperature gradient.

Sequestration (Pumping) Rate
The expected concentration of subsurface H 2 O may be lower than the equilibrium concentration, due to the slow migration at low temperatures. Here the sequestration rate is estimated with a random-walk model.
Diffusion and adsorption are modeled in a one-dimensional setting where molecules move up or down in discrete steps and a soil grain is reduced to a site on a discrete vertical grid. Every incoming molecule is assigned a residence time τ, picked from a probability distribution such that 1/τ is distributed exponentially and where the time average is over a model time interval Δt. The molecule then leaves the grain surface at thermal speed. In the temperature range explored, the flight time for the short intergrain distances is negligible compared to the residence times.
Temperatures are advanced with time steps Δt, chosen as 1/96 of a lunation. Within each such time interval, molecules can jump repeatedly until the next departure time is beyond t n+1 = t n + Δt. The adsorbate concentration θ is also updated in intervals of Δt. When more than one monolayer's worth of computational particles is present, molecules that are covered by other molecules are considered immobile.
The model calculations are initialized with a completely dry domain. At each time step, one computational particle is added on the surface. Those that move upward from the surface are assumed lost. A delivery rate of 1000 kg m −2 Gyr −1 corresponds to about 3 monolayers per year or 1/4 of a monolayer per month. Hence, one monolayer is represented by 96 × 4 = 384 computational particles. As time progresses, the model keeps track of the locations and residence times of more and more particles. The mean free path is set to λ = 100 μm, which is of the same order of magnitude as the grain size of lunar dust in Apollo samples (McKay et al. 1991). Figure 4 shows the results for the amount sequestered at the site with the lowest mean temperature among the three study sites. The development of a hydrated layer separated from the surface by an adsorbate-poor layer is apparent (Figure 4(a)). The maximum adsorbate concentration first exceeds a monolayer after 23 kyr. The pumping pushes molecules into the slowly hydrating layer, where the molecules are barely mobile. The depth to the maximum H 2 O density is expected to slowly increase with time toward the asymptotic depth for ice.
The total amount sequestered (Figure 4(b)) increases over time at a rate of about 1.2 μg m −2 yr −1 . As a fraction of the amount delivered to the surface (Figure 4(c)), the columnintegrated sequestered mass is 0.13% over the first 100 kyr but decreases with time, as more and more molecules migrate back up from the growing ice reservoir. Extrapolating to 1 Gyr, the fraction may be 0.05% or 0.5 kg m −2 Gyr −1 . The Moon has been in its current spin axis orientation probably for the past 2-3 Gyr (Arnold 1979). Overall, a column-integrated abundance on the order of 1 kg m −2 (the equivalent of a 1 mm thick layer) is expected for the aforementioned rate of supply to the surface. If the amount of water delivered to the polar regions is smaller than assumed here, the abundance in the subsurface will also be lower. Schorghofer & Aharonson (2014) estimated that the amount sequestered is λ/Δz times the supply rate, where Δz is the depth to where ice starts to accumulate, which in this example is roughly 10 cm. This yields 0.1%, the correct order of magnitude. Reiss et al. (2021) calculated the amount sequestered over the first 300 lunations, using a lower supply rate. At their study site "P2," 0.05% of the delivered amount is sequestered. This is remarkably similar to the result obtained here, given that the site, as well as many aspects of the models, are different.
After 100 kyr, a centimeter-thick icy layer has formed (Figure 4(a)). Whether the thickness of the hydrated layer reaches the full few decimeters of the equilibrium solution (Figure 3(d)) is unclear. In addition, impact mixing will interfere with the subsurface migration and deplete some of the H 2 O while protecting another portion. Impact mixing is estimated to turn over about 10 cm Gyr −1 (Gault et al. 1974). Roughly estimated, 1 kg m −2 of H 2 O distributed over about 10 cm depth represents an average density of 10 kg m −3 , nearly 1% by mass of particulate silicate rock.
The rate of sequestration will also depend on grain properties. Larger grains have a smaller specific surface area but also larger pore spaces and therefore a longer mean free path. Assessments of the grain size dependence and the role of surface microtopography (Davidsson & Hosseini 2021;Sarantos & Tsavachidis 2021) are postponed to a later study.

Discussion: Relation to Observations
A predicted concentration on the order of 1 wt% far exceeds the trace amounts of nonvolatile hydrogenous species found in lunar soils, as measured in Apollo samples (McCubbin et al. 2015) or through remote sensing (Honniball et al. 2021), which are typically on the order of 1-1000 ppm. The amount could be easily detected with upcoming neutron spectrometers on rovers, such as those carried by the Volatiles Investigating Polar Exploration Rover and MoonRanger (Schweitzer et al. 2021).
Neutron spectroscopy from lunar orbit has revealed a hydrogen excess in both polar regions, and the different response of thermal, epithermal, and fast neutrons suggests the hydrated layer is buried (Feldman et al. 1998;Lawrence et al. 2006), perhaps by 10 cm. This is consistent with the predicted thickness of the desiccated layer (Figures 3(d) and 4(a)). Boynton et al. (2012) divided the H population into localized neutron suppression regions and a latitude-dependent H enhancement. This second population could be due to the subsurface hydrated layer. The areas where vapor pumping occurs is more than five times larger than the area of cold traps (Schorghofer & Williams 2020) and therefore has a broad geographic distribution. Gläser et al. (2021) identified areas in the north polar region with buried excess hydrogen at temperature conditions favorable for vapor pumping.
Measurements by the Lunar Atmosphere Dust and Environment Explorer (LADEE) revealed that the lunar surface releases a water-group species (OH or H 2 O) during meteoroid streams but not from the more time-continuous bombardment by dust-size particles (Benna et al. 2019). Combining LADEE observations with modeling suggests that the minimum mass of impactors needs to be about 0.15 g (Benna et al. 2019), which corresponds to a projectile of about 0.5 cm in diameter. Excavation or disturbance takes place to a depth multiple times the size of the impactor. Hence, one can hypothesize that centimeter-size impactors disturb the hydrated layer, but millimeter-size impactors do not. Subsurface hydrated layers are restricted to the polar regions, but the released water would quickly spread globally through exospheric transport and could be observed even at the equatorial latitude of the LADEE spacecraft.

Conclusions
The sequestration mechanism requires two components: A sufficiently cold subsurface and temperature cycles that accelerate diffusion into the subsurface.
If a 1 nm thick layer of H 2 O is delivered to the surface each year, then a porous subsurface layer continuously colder than about 114 K retains molecules indefinitely, resulting in the buildup of ice in the subsurface. Landed missions exploring for ice outside of cold traps should target these temperature conditions. At higher subsurface temperatures, only adsorption is expected to occur, with an abundance that quickly diminishes with the difference from the threshold temperature. This threshold behavior can provide a measurement of the long-term delivery rate of H 2 O.
A desiccated layer overlies the hydrated layer. The thickness of the desiccated layer is set by the thermal skin depth, which in lunar dust is typically centimeters. This potentially provides an explanation for why centimeter-size meteoroid impacts release water-group species from the lunar surface, whereas millimeter-size meteoroids do not (Benna et al. 2019). A desiccated layer is also inferred from neutron counts (e.g., Feldman et al. 1998;Lawrence et al. 2006).
Software: Source code for the models is available online (Schorghofer 2021).