Mapping of Ice Storage Processes on the Moon with Time-dependent Temperatures

Lunar cold traps are defined by extremely low time-integrated sublimation loss, so that water ice is expected to accumulate in them. Due to the strong dependence of the sublimation rate on temperature, they have heretofore been delineated by the peak rather than the time-averaged sublimation rate. Here, time-averaged sublimation rates are calculated for the south polar region of the Moon based on 11 yr of Diviner surface temperature measurements. Data for each spatial pixel are binned according to subsolar (diurnal) and ecliptic (seasonal) longitude. Frequency-domain filtering of temperature time series is applied to interpolate and smooth the data set. The cold trap area is 17,000 km2 from 80°S to the south pole, about 24% larger than defined by peak temperature. In addition, subsurface ice stability is mapped based on solutions of the heat equation with Diviner surface temperatures as the upper boundary condition. Even a thin layer of dust reduces the sublimation loss dramatically. A third potential mechanism for ice storage, vapor pumping by temperature cycles, is also mapped, based on a model for the time-variable population of adsorbed water molecules on the lunar surface. The area of significant pumping is estimated to be 96,000 km2 poleward of 80°S. In summary, a new processing technique has been developed that exploits information about the time variability of lunar surface temperatures measured by Diviner, and it results in new maps for the geographic distribution of potential water-ice deposits in the south polar region of the Moon.


Introduction
The distribution and abundance of water ice in the lunar polar regions remain some of the most important and unanswered questions in lunar science (e.g., Hayne et al. 2014;Elvis et al. 2016). Classical theory suggests ice should reside within cold traps (Watson et al. 1961a(Watson et al. , 1961bArnold 1979), but the observational evidence points toward a more complex situation (Hayne et al. 2014). The retention of volatiles, including H 2 O, depends sensitively on surface or near-surface temperature. The extent of cold traps is commonly defined by a peak temperature threshold (110 K), whereas it is fundamentally defined by the sublimation rate (0.1 m Gyr −1 ). Here, lunar cold traps are mapped based on time-integrated sublimation rates instead of peak temperature. Diviner bolometric temperatures over 11 yr are processed to calculate time-dependent quantities. This approach is then extended to more accurately map the stability of buried ice and "thermally pumped" ice, two proposed alternative ice storage mechanisms (Paige et al. 2010a;Schorghofer & Aharonson 2014).
Cold traps lie within permanently shadowed regions (PSRs), but not all PSR areas act as cold traps, because terrain irradiance contributes to the surface energy balance (Hayne et al. 2020). The Diviner Lunar Radiometer Experiment, on the Lunar Reconnaissance Orbiter (LRO), measures surface temperature at high spatial resolution (Paige et al. 2010a(Paige et al. , 2010b. Previously, maps of lunar cold traps based on Diviner temperatures were published in Paige et al. (2010a), Hayne et al. (2015), and Williams et al. (2019). Williams et al. (2019) identified the draconic period in the surface temperatures and produced winter and summer maps of cold traps. The draconic year is the time between successive peaks in the subsolar lunar latitude. It is shorter than Earth's year because of the precession of the Moon's orbit around Earth. Over the draconic year (346.62 Earth days), the decl. of the Sun as seen from the Moon (the subsolar latitude) changes between ±1°.57. In addition to the synodic month (29.53 days) and the draconic year, the Moon-Sun distance and the 18.6yr precession cycle influence the solar energy input, but the amplitude of these cycles is small and not apparent in the Diviner data set.
In terrestrial geography, the "potential evaporation" (PE) or "potential evapotranspiration" (PET) is the evaporation rate if water is present (Mather 2009). For the Moon, we can define a "potential sublimation rate." The abbreviation PSR is already reserved for "permanently shadowed region," so we will refer to the (long-term average of the) potential sublimation rate as PE as well. In terrestrial hydrology, it is common to compare the annual potential evaporation with the annual precipitation (P). The "aridity index" is defined as P/PET. Hyperarid zones (extreme deserts) are defined by P/PET < 0.05 (Mather 2009). Lunar cold traps are areas where the influx of water molecules is larger than PE. The influx rate P is not well constrained for the Moon; a plausible value is 1m Gyr −1 (Arnold 1979). But the potential sublimation rate, PE, depends very strongly on temperature, such that the geographic extent of cold traps is only weakly dependent on P (or on the chosen threshold for PE), which justifies defining cold traps based purely on a PE threshold.
The sublimation rate of ice in vacuum is given by the Hertz-Knudsen formula (Kramers & Stemerding 1951;Persad & Ward 2016) where p s is the equilibrium vapor pressure of ice, μ the molar mass (18.015), R u the universal gas constant, and T the temperature. The main temperature dependence lies in p s (T). The sublimation coefficient, 0α1, is close to unity at low temperature (Haynes et al. 1992). Criteria for defining cold traps are The first two criteria, Equations (2) and (3), are equivalent and correspond to the classical approach of using the peak temperature to define cold traps. Criterion (4) is based on the time-averaged potential sublimation rate and will be used in this work. Cold trapping theory is consistent with observations of ice on Mercury (Neumann et al. 2013;Paige et al. 2013;Chabot et al. 2014) and on Ceres (Ermakov et al. 2017). If classical cold trapping also operates on the Moon, the cold traps determine the distribution of ice in the lunar polar regions. Proxy measurements of H 2 O ice content in the near-surface of the Moon provide conflicting results for the geographic distribution of ice, but some of the measurements are consistent with all cold traps having stored ice (e.g., Feldman et al. 2001;Colaprete et al. 2010;Rubanenko et al. 2019). Other processes also influence the presence of ice on the Moon, such as meteoric impact vaporization, sputtering, and meteoric impact ejection (e.g., Farrell et al. 2019;Costello et al. 2020), with the unknown supply rate of H 2 O being perhaps the greatest uncertainty. Nevertheless, temperature is an important (and arguably the most important) factor that governs ice storage on the Moon.

Diviner Data Reduction and Binning
The Diviner Lunar Radiometer Experiment is an infrared and solar radiometer aboard NASA's Lunar Reconnaissance Orbiter (LRO; Paige et al. 2010b) and has continuously operated since July 2009. It has a spectral range of 0.35-400μm in nine spectral channels. Each channel consists of a line array of 21 detectors nominally operating in a nadirpointing push broom configuration with an integration period of 128ms. Variations in LRO's altitude have been minimal over the south polar region, and as a result, Diviner's surface footprint has remained relatively fixed in the region during the mission with a cross-track width of ∼140m and an in-track length of ∼400 m.
As part of this study, Diviner Reduced Data Records (RDR) data from 2009 July 5 to 2020 March 15 was compiled into seasonal and fixed subsolar longitude bins. The calibrated radiance measurements from Diviner channels 3-9 covering wavelengths from 7.55 to 400μm were map-projected onto a polar stereographic grid from 80°S to 90°S at a resolution of 240 m/pixel and converted to bolometric brightness temperatures (Paige et al. 2010a). Observations with emission angles >10°, or acquired during off-nadir activities, were excluded as a reduction in emitted radiance with increasing emission angles has been observed (Bandfield et al. 2015;Warren et al. 2019). Prior to binning, the RDR measurement locations were updated to account for topography and the individual effective field of view of the detectors, which were modeled for each RDR record to account for the spacecraft motion and the thermal response time of the detectors as described in Williams et al. (2016) and Sefton-Nash et al. (2017).
The radiances are then partitioned into bins of subsolar longitude (λ s ) and ecliptic longitude (L s ) to characterize the diurnal and seasonal temperatures, where L s is defined by the R.A. of the apparent position of the Sun as seen from the Moon in the reference frame of the ecliptic north pole and the Moon's pole. The L s calculations were done using the Jet Propulsion Laboratory's Navigation and Ancillary Information Facility (NAIF) WebGeocalc tool (Acton et al. 2017). Defined this way, L s provides the angular description (0°-360°) of the draconic time of year, with 0°being when the Moon crosses the ascending node (i.e., spring equinox in the northern hemisphere).
A temperature map is produced for each diurnal and seasonal bin. Williams et al. (2019) have created 96×2 maps. If the number of bins is too large, most bins will contain no data. If they are too few, information about the time dependence of the temperature is lost. A relevant question is how to partition the available data among monthly and seasonal bins. The Appendix considers this choice for an idealized temperature dependence. The optimal number of bins turns out to be proportional to the square root of the respective temperature amplitude. The Appendix also provides an estimate of the ratio of monthly and seasonal amplitudes, which is about 6 in the polar region. In some areas, the influence of the season is larger than that of the lunation. This analysis motivated us to move from 96×2 bins to a smaller number of diurnal bins in favor of a larger number of seasonal bins.
Our maps are produced based on data organized in 24×6 time bins, where the first seasonal bin stretches from L s =0°-60°. The second seasonal bin is centered on L s =90°, when the Sun reaches its maximum decl. of 1°.57 and is at its lowest elevation in the southern hemisphere. With this choice for bin centers, the first three bins correspond to southern winter (decl. >0°) and the other three bins to southern summer (decl.<0°). Diurnal bins are based on the subsolar longitude λ s . The relation to local time (in hours) is where λ is the geographic longitude in degrees (positive east). Figure 1 shows the spatial distribution of the time coverage. The median of the time coverage is 91 out of 144 bins (63%). In the 96×2 data set by Williams et al. (2019), the coverage was 121 out of 192 (also 63%). No pixel has all of its 144 time bins populated with measurements.
Maps are in polar stereographic projection from 80°S to 90°S , based on a radius of R=1737.4 km, which results in maps with 2533×2533 pixels. The area of a pixel in a polar stereographic projection is calculated as where 1/k 0 =240 m and x and y are the pixel coordinates relative to the center pixel.

Two-dimensional Time-domain Interpolation (Gap Filling)
Whereas a maximum temperature can be determined even when many time bins are empty, calculations of time-integrated sublimation rates require that all bins are populated. To fill empty bins, the 24×6 time bins are downsampled individually for each pixel rather than globally for the entire map. Figure 2(a) shows the temperatures in the 24×6 time bins for one pixel. About one-third of the bins are empty (which is typical), and they are filled with the downsampled version of the same data. The same temperature data with downsampled bins are shown in panel (b) of the figure. From 24×6 bins, versions with 12×2 bins (2 × 3 boxes) and, if necessary, 6×1 bins (4 × 6 boxes) are generated. Panel (c) shows the original and filled data combined, without any remaining data gaps. If no occupied bin is found within a 4×6 downsampled time box, no sublimation rates are calculated for this pixel, but this only occurs near the margin of the geographic domain.
Figures 2(e)-(h) shows another example of the bin filling procedure, this time from within a PSR where the seasonal amplitude is clearly discernible. The three seasonal bins averaged at the downsampling step (Figures 2(b), (f)) correspond to the winter and summer seasons (as defined by solar decl.).
Reasonable interpolations are obtained even when more than half of the bins are missing. The missing bins are not concentrated at a particular local time or a particular season (Figures 2(a), (e)), and populated bins are nearby (in the time domain), which makes this approach to interpolation effective. When fewer than one-third of the bins are populated with original data, no values are calculated for the pixel, but this only occurs on the very margin of the domain (around 80°S).

Frequency-domain Filtering
Missing data can be interpolated and existing data smoothed by filtering in the Fourier domain. The Fourier domain is a natural choice, because the signal is periodic. To be able to perform a Fast Fourier Transform of the time series, missing values are first replaced with interpolated values, as described above (Section 2.2). Then, the Fourier transform is performed along the diurnal time direction, a low-pass filter is applied, and the inverse Fourier transform is taken.
We experimented with various filter shapes and settled on the Welch window. With N time bins, the Fourier frequencies f range from 0 to N/2 (with the time period normalized to 2π). A Welch (parabolic) window is of the form = - ) . Examples of the results of the frequency-domain filtering are shown in Figure 3. Full symbols represent bins populated with temperature measurements, and empty symbols are those created after two-dimensional time interpolation and onedimensional Fourier filtering. For example, in Figure 3(c), erratic up-and-downs in the data are smoother after Fourier filtering. We find the Fourier-filtered temperatures to be so robust that we use them not only to populate empty bins but also to replace measured temperatures, which reduces the scatter in the data. Below we refer to this as the "smoothed" version of the data set, as opposed to the unsmoothed version where time bins originally populated with temperature measurements are never overwritten.
When the bin with the (presumably) maximum temperature is missing any measurements, the interpolation scheme never produces a new maximum, but the Fourier filtering can produce a temperature value larger than the temperature values in all bins that are populated with the measured data.

Sublimation Rates
Murphy & Koop (2005) reviewed vapor pressure data for crystalline (hexagonal) ice, which is the crystal structure of all natural snow and ice on Earth. All of the commonly used expressions for the vapor pressure of ice, except that of Marti & Mauersberger (1993), are within 1% of each other for temperatures between 170 and 273K, including that derived from the Clapeyron equation with a temperature-independent sublimation energy. Here, we adopt Equation (7) ( ) The sublimation rate of ice in vacuum is given by Equation (1). The coefficients in Equation (7)   A fit to the inverse of Equation (8) is The temperature is then refined with a few iterations of the Newton method, using the analytical derivative of Equation (8).
As above, E is here in units of kg m −2 s −1 and T in units of Kelvin. (We convert seconds to Gyr using a year length of 365.24 days.) At 110K, the sublimation rate of ice is too low to be measurable in the laboratory, but the saturation vapor pressure can be extrapolated to lower temperatures based on the latent heat. Theoretical expressions for crystalline ice in the form of Equation (7) can be robustly extrapolated to lower temperatures (Jancso et al. 1970;Murphy & Koop 2005).
However, amorphous ice starts to form below about 140K (Sack & Baragiola 1993). The transition temperature depends on the rate of deposition and is lower at lower deposition rates. Low-temperature ice can exist in three amorphous and two crystalline forms (Jenniskens et al. 1998). Hence, the sublimation rate of low-temperature ice depends on more than temperature, resulting in an intrinsic uncertainty about the pertinent sublimation rate. Adding to the uncertainty is that ice could have been deposited in lunar cold traps at a different temperature from their current environment, so it could be crystalline or amorphous in nature. The vapor pressure of amorphous ice is not necessarily much higher than that of crystalline ice. For example, Speedy et al. (1996) found roughly a factor of 2 difference in desorption rates, but Kouchi (1987) measured a vapor pressure increase by an order of Figure 2. Illustration of the data gap filling procedure. The data in each spatial pixel is divided into 24 diurnal bins and 6 seasonal bins (panel a, where black indicates missing bins). A downsampled version of the data is created with coarser bins (panel b), and these values are used to fill missing bins (panel c). Panels (a)-(d) are for a pixel at 81.43°S 45°W. Panels (e)-(h) correspond to a pixel within the PSR of the Haworth crater (86°. 90S 3°. 95W). Diurnal bins 1-24 are based on subsolar longitude. Seasonal bins 1-3 stretch from L s =0°-180°(southern winter) and bins 4-6 from L s =180°-360°(southern summer). magnitude relative to crystalline ice. For comparison, the ratio of sublimation rates at 120K and 110K, as extrapolated for crystalline hexagonal ice, is a factor of nearly 100. Hence, a pessimistic scenario of an order of magnitude uncertainty in the sublimation rate corresponds to a 5K difference in temperature. Table 1 shows commonly chosen threshold values with corresponding values for temperature, sublimation rate, and molecular residence time. The latter is calculated as q m E 0 , where q = 10 0 19 m −2 is the number of H 2 O molecules per monolayer, which follows from the density of ice. This is the residence time for a H 2 O molecule on ice; adsorbed molecules potentially have much longer residence times.
The instantaneous surface temperature can be above the triple point temperature of H 2 O. These locations are far from suitable for ice stability but are evaluated to create complete maps. At temperatures above 273.16K, the potential evaporation rate is calculated based on the vapor pressure of liquid water (Hyland & Wexler 1983;Murphy & Koop 2005).

Figures 4(A) and (B)
show the maximum potential sublimation rate, E max t ( ), and the time-averaged potential sublimation rate, mean t (E), in the south polar region. They vary by more than 40 orders of magnitude over the area shown, and the color scale is truncated at both ends of the range. Figure 4(C) shows the ratio of max t (E) to mean t (E). Typically, E max t ( ) is larger than mean t (E) by an order of magnitude. The change in cold trap contours (for 100 kg m −2 Gyr −1 or 114 K) is displayed in Figure 4(D) and shows that cold traps are larger and more numerous for mean t (E). Figure 5 shows cold trap areas for various thresholds and methods. It compares results for max t (E) with mean t (E) as well as the smoothed and unsmoothed data sets. For the smoothed data sets, measured temperatures are replaced with their Fourier-filtered values, whereas in the unsmoothed data set, only unpopulated time bins are replaced with the Fourierfiltered values. Using a threshold temperature of 110K (114 kg m −2 Gyr −1 ), the cold trap area as determined by the unsmoothed maximum temperature is 13.1×10 3 km 2 , consistent with the 13.0×10 3 km 2 determined by Williams et al. (2019). Using the smoothed maximum temperature results in a slightly larger area of 14.1×10 3 km 2 . Using mean t (E) instead of max t (E) results in significantly larger areas of 16.4×10 3 km 2 , unsmoothed, and 17.5×10 3 km 2 , smoothed. Because the physically appropriate definition of a cold trap is through mean t (E), not max t (E), the larger area is more appropriate. For the smoothed and unsmoothed data sets, the cold trap area based on mean t (E) is 24% larger than that based on E max t ( ). Changing the evaluation criterion from max t (E) to mean t (E) makes more of a difference to the cold trap area than changing E thres from 10 to 1000kg m −2 Gyr −1 .
Statistical fluctuations of temperature values have the tendency to increase the peak temperature, as a statistical deviation below the true value may be superseded by the second highest temperature. Hence, because of outliers, the Fourier-smoothed version is expected to be more accurate than the unsmoothed version. Overall, 17×10 3 km 2 is our best estimate for the cold trap area poleward of 80°S when delineated by a sublimation rate of 114kg m −2 Gyr −1 and, within two significant digits, also for 100kg m −2 Gyr −1 .
Poleward of 80°S, there are 32 cold traps larger than 100 km 2 , 171 larger than 10km 2 , and 703 larger than 1km 2 , based on a threshold of 100kg m −2 Gyr −1 and the smoothed data set. Using the peak sublimation rate instead, these numbers would be 18, 129, and 562, respectively. The largest cold traps are the floor of the Shoemaker crater (1314 km 2 ), the floor of the Haworth crater (1297 km 2 ), and the northern portion of the Amundsen crater (1186 km 2 ).
The aforementioned aridity index, defined as P/PE= E thres /mean t (E), is shown in Figure 6 for a supply rate of P=1000kg m −2 Gyr −1 . This mirrors the map of PE in Figure 4(B) but emphasizes the extreme desiccation in most areas and the contrast to water-accumulating areas (the cold traps).

Subsurface Ice Stability
For subsurface ice, the sublimation loss is reduced because the overlying ice-free layer acts as a diffusion barrier, but also because the ice does not experience the same peak temperature as the surface. The loss rate of subsurface ice, E ss , in terms of the loss rate of ice exposed on the surface, E, is roughly (Schorghofer & Taylor 2007;Schorghofer 2008 Hayne et al. (2017), which corresponds to the "Diviner standard thermal model." Specifically, we use the same values and depth-dependence for thermal conductivity, heat capacity, and regolith density, and a gradual transition of thermal properties with a characteristic depth of H=7cm. The geothermal heat flux at the bottom of the computational domain is set to 0.018Wm −2 (Langseth et al. 1976). The thermal properties are considered to be time independent, and the specific heat capacity is evaluated at the mean surface temperature. The temperature profile is initialized with the temperature profile expected based on the mean surface temperature, thermal conductivity profile, and geothermal heat flux.
The input surface temperature time series has two periods (the synodic month and the draconic year). The ratio of their periods is about 11.74, so the periods are incommensurate (not a ratio of small integers) and the time series does not repeat for a long time. For every synodic month, 24 time steps are taken.
After the passage of one draconic year, surface temperature values for the next L s bin are used. The total run stretches 70 months. Time-averaged sublimation rates are calculated at each depth from the last 12 synodic months, which capture slightly more than one draconic year. The thermal skin depth characterizes the penetration depth of thermal cycles but is properly defined only with uniform The nominal skin depth of the synodic month is 5cm and the nominal skin depth of the draconic year is 19cm. Large cohesive rocks increase thermal conduction, and the thermal wave will penetrate deeper.
Due to the decaying temperature amplitude with depth, the mean t (E(z)) decreases with depth, but due to the geothermal heat, it will increase beyond the influence of the thermal wave. The shallowest depth with mean < E E t ss thres ( ) is determined, if there is any. The depth of the computational domain was chosen to be 1m.
Vertical temperature profiles for one location are shown in Figure 7. At this location, ice is not stable at any depth. In the top decimeter, the mean E t ss ( ) changes by nine orders of magnitude. Even a thin overlying layer dramatically reduces the loss rate. This is in part due to the attenuation of the peak temperature and in part due to the diffusion barrier.

Results for Subsurface Ice Stability
Figure 8(A) shows the mean surface temperature, where the time average is obtained by weighing each time bin equally. In other words, it is compensated for nonuniform time coverage. The temperature averaged over this area is 113K, so subsurface ice stability can be expected to be extensive. As mentioned above, at 1m depth, the threshold temperature increases to around 131K, and the mean surface temperature approximately equals the temperature at 1m depth. Figure 8(B) shows the estimated depth to 100kg m −2 Gyr −1 loss. Subsurface ice is stable in most of the region, as was already discovered by Paige et al. (2010a). Based on the subsurface profile shown in Figure 7, such extensive and shallow depths are physically reasonable. However, ice delivered from above would have had to be buried quickly; otherwise, it experienced significant loss before it was buried. These calculations only determine the stability of ice, if ice preexists at this depth. Figure 9 shows the area where subsurface ice is stable as a function of depth. The area initially increases quickly with depth, mostly due to the attenuation of the peak temperature.
The area of stability doubles at 2cm depth. If ice is not stable in the top 20cm, then it is rarely stable at greater depth. The area considered is limited to the region poleward of 80°S and would be larger if more equatorward regions were added. The results shown in Figure 9 also illustrate that the area with stable subsurface ice is not sensitive to the choice of the molecular mean free path ℓ or E thres .

Vapor Pumping Calculation Method
Temperature cycles can drive volatile water molecules on the surface into the regolith. This "pumping" effect is expected to operate on the lunar surface, typically on pole-facing slopes near to but outside of cold traps, although at low yield (Schorghofer & Aharonson 2014). Quantifying this effect requires a model for the population of volatile H 2 O molecules on the surface, as well as subsurface temperatures. The surface population model is exactly the same as in Schorghofer & Aharonson (2014). Each adsorbed molecule has a residence time that depends on temperature and the areal density of water molecules. The flux of molecules leaving the surface is is the sublimation rate of pure ice, and 0f1. The function f is based on adsorption measurements on a lunar sample at a single temperature (Cadenhead & Stetter 1974). The function E is for amorphous ice according to Sack & Baragiola (1993), slightly different from the parameterization for crystalline ice used above. The population density θ is then determined by the differential equation The supply, s, is assumed to be at a uniform rate of 10 12 molecules m -2 s, which according to Table 1 corresponds  to 1m Gyr −1 . The destruction rate by space weathering, w, is usually negligible but is included. The ordinary differential Equation (12) involves terms that change at different timescales and is solved numerically with an integrator for stiff equations, using a backward time difference instead of a forward time difference. Even so, a very small time step is necessary to achieve mass conservation. We chose 12,000 time steps per lunation (213 s). The temperature time series is generated as described in Section 4.1, covers 70months, and time averages are taken over the last 12 synodic months. The pumping differential is defined by Net accumulation occurs when ΔE>0 and the mean t is the long-term average. The last term is the sublimation rate at the shallowest depth with ice, and for the purpose of determining whether or not ice accumulates at any depth, the minimum over depth can be taken. For simplicity, we will approximate » E E T min mean mean z t t ( ( )) ( ( )). The pumping differential ΔE has the same units as the sublimation rate, but the amount that accumulates in the subsurface through pumping is much less than ΔE (Schorghofer & Aharonson 2014). Here, we only map ΔE, which is bounded from above by the supply rate s.

Results of Pumping Calculations
Pumping crucially relies on the time dependence of temperatures and does not occur in a static thermal environment. Volatile H 2 O molecules on the surface can be lost to space but a fraction will diffuse into the regolith. A map of the pumping differential has been published in Schorghofer & Aharonson (2014), using a time dependence parameterized by two values, the mean and maximum temperature. The time dependence used here is far more detailed and parameterized by typically 91 independent parameters (Section 2.1). This is also the first time that pumping by seasonal cycles, which penetrate more deeply than the monthly cycles, is evaluated.
A map of positive pumping differentials ΔE is shown in Figure 10 and resembles the one published in Schorghofer & Aharonson (2014). Pumping is strongest in areas adjacent to cold traps. Within cold traps, pumping is positive, but very weak. The total area where any pumping occurs is 121×10 3 km 2 poleward of 80°S. The pumping is substantial (ΔE>0.5s) over an area of 96×10 3 km 2 , more than five times the area of surface cold traps.
These model results involve significant uncertainties due to uncertainty in the input parameters, such as H 2 O infall rates, adsorption residence times, and the role of surface roughness. The amount of subsurface ice that accumulated is proportional to ΔE, the molecular mean free path, and time.

Summary
The Diviner data set provides lunar surface temperature measurements at high spatial resolution and with extensive time coverage. A data processing method has been developed that generates quasi-continuous time series from this large data set. Temperature measurements from each spatial location were binned according to the subsolar longitude (related to local time) and ecliptic longitude (season). Empty time bins are filled by two-dimensional time-domain interpolation using box averages. To further improve the interpolation and reduce noise in the data, frequency-domain filtering is applied.
The map of south polar cold traps, properly defined by the long-term potential sublimation rate, is shown in Figures 4(B) and 4(D), and 6. The total cold trap area is larger by onequarter compared to the same data set evaluated with the traditional peak temperature criterion. Given that all cold trap   areas may be ice-bearing, this adds to the regions of interest for the future exploration of the south polar region.
Subsurface thermal conditions were determined using the same time-continuous model of Diviner temperature cycles combined with numerical solutions for subsurface heat conduction. They capture the effect of monthly and seasonal temperature variations as well as that of a diffusion barrier. A protective layer of regolith reduces sublimation loss dramatically, and areas where subsurface ice would be stable are ubiquitous in the south polar region (Figure 8(B)). The area of ice stability doubles at 2cm depth, but ice would have to be buried quickly to survive. The model calculations also reveal that if ice is stable (against sublimation loss) at any depth, it usually is already stable within 20cm depth, the typical penetration depth of the draconic temperature cycle. Pumping of water molecules into the porous regolith occurs over areas much larger than the surface cold traps (Figure 10).
Averaged over the offset, the error due to the time variation of temperature, in the bin with the maximum temperature, A 7 2 ( ) ( ) Figure 11. Illustration of the time variation within a single bin close to the maximum. Here the bins are 1 local hour wide.
The error is negative, i.e., it underestimates the maximum, as evident from Figure 11. This result can be applied to the lunation, with ecliptic longitude fixed, or to the draconic year, with local time or subsolar longitude fixed. The average of temperature values within a bin is lower than the true maximum. In terms of the number of bins, For A=100K and B=24, this would only be −0.6K. With only two bins, however, the deviation is significant.

A.2. Optimal Number of Bins with Two Periods
The total error is The first term is the statistical error, and the other two terms arise from the time variability within a bin, Equation ( The relative number of bins is determined solely by the ratio of temperature amplitudes. At the equator of the Moon A L ? A D , but in the polar regions, the seasonal temperature amplitude A D can be large.

A.3. Estimates of Temperature Amplitudes
Temperature amplitudes in the south polar region are estimated based on the hourly temperature maps with six seasonal bins (24 × 6). The monthly max-min temperature difference is calculated as the mean max-min difference among the six seasonal bins, and the seasonal max-min difference is calculated as the mean difference along the 24 diurnal bins. Only pixels where at least half the bins are populated are considered. The amplitude ratios are shown in Figure 12.
The median of the ratios determined by this method is 6.2 in the south polar region. In the vicinity of the pole, the ratio tends to be below this median value. In some areas, the seasonal amplitude is even larger than the monthly amplitude (e.g., in the Haworth crater). With » 24 6.2 10, a desirable choice would be 24×10 bins, but these are more bins than Figure 12. Estimates of the ratio of diurnal to seasonal temperature amplitude in the south polar region. The ratio is evaluated only for pixels with ample time coverage, and the color map is logarithmic.