Lunar North Polar Cold Traps Based on Diurnally and Seasonally Varying Temperatures

Lunar cold traps are defined by extremely low sublimation rates, such that water ice could have accumulated in them. Here time-averaged sublimation rates are calculated for the north polar region of the Moon based on over 14 years of Diviner surface temperature measurements. Data for each spatial pixel are binned according to subsolar (diurnal) and ecliptic (seasonal) longitude. The cold trap area poleward of 80°N is about 32% larger when defined by a time-average sublimation rate instead of by peak temperature. Apparently sunlit cold traps are identified, e.g., in Lenard Crater, where modeling of direct illumination reveals that the Sun briefly rises above the horizon each Draconic year. The true cold trap area is smaller than what can be determined from Diviner data. Also presented are north polar maps for the potential sublimation rate of relic buried ice and for subsurface cold trapping.


Introduction
Water ice is thought to have accumulated in the lunar polar regions in cold traps, where temperatures are so low that the sublimation rate of ice into vacuum becomes negligible (Arnold 1979;Lawrence 2017).The Diviner Lunar Radiometer Experiment on board the Lunar Reconnaissance Orbiter (LRO) has provided brightness temperatures of the lunar surface (Paige et al. 2010a;Williams et al. 2017) and resulted in detailed maps of lunar cold traps based on peak temperatures (Paige et al. 2010b;Hayne et al. 2015;Williams et al. 2019;Landis et al. 2022).Instead of using peak temperature, Schorghofer & Williams (2020) produced cold trap maps based on diurnally and seasonally varying sublimation rates for the south polar region.Here, lunar cold traps are mapped in the same way for the north polar region (Section 3).
Cold traps usually lie within permanently shadowed regions (PSRs), because direct sunlight is nearly always sufficient to raise the surface temperature above the stability threshold.The total cold trap area is known to be smaller than the total PSR area (Hayne et al. 2021), because terrain irradiance contributes to the surface energy balance.However, this is not necessarily what is observed from Diviner data at every location.Comparison of the extent of individual PSRs with individual cold traps leads to the identification of apparently sunlit cold traps, but these are explained by undersampling (Section 4).
The same set of Diviner temperature maps is also used to calculate subsurface temperatures and time-averaged subsurface sublimation rates, leading to new maps of two types of subsurface stability in the north polar region (Section 5).

Diviner Data Reduction and Processing
The Diviner Lunar Radiometer Experiment is an infrared and solar radiometer on board NASA's Lunar Reconnaissance Orbiter with a spectral range of 0.35-400 μm in nine spectral channels (Paige et al. 2010a).Diviner has operated nearcontinuously since 2009 July.Interruptions to the instrument's nominal nadir mapping have occurred for routine calibration sequences and intermittent off-nadir observations and other special spacecraft operations and activities.
After an initial ∼10 week commissioning phase in an elliptical orbit, LRO entered a ∼50 km near-circular mapping orbit followed by a transition back to an elliptical orbit at the end of 2011 with an apoapsis over the north polar region at an altitude >100 km.The inclination of the spacecraft was initially ∼90°providing complete spatial coverage of the polar regions.However, the inclination of the orbit has migrated to ∼85°latitude during the mission lifetime due to the Moon's 18.6 yr nodal precession.This has generated an expanding coverage gap around the poles while increasing the density of ground tracks over lower-latitude polar regions (Petro et al. 2022).
As part of this study, Diviner Reduced Data Records (RDR) data from 2009 July 5 to 2024 January 20 (more than 15 Draconic years) were 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°N to 90°N at a resolution of 240 m px -1 and converted to bolometric brightness temperatures (Paige et al. 2010b).Observations with emission angles >10°, or acquired during off-nadir activities, were excluded to minimize phase angle effects that have been observed to result in a reduction in emitted radiance with increasing emission angles (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, Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 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 Draconic year (346.62Earth days) is the time between successive peaks in the subsolar lunar latitude.(The lunar pole points toward the constellation Draco.)Radiances are partitioned into bins of subsolar longitude (λ s ) and subsolar latitude to characterize the diurnal and seasonal temperatures.Over long periods, the average time spent in each bin is the same (due to the 18.6 yr precession cycle), so the long-term time average is the arithmetic average over these temporal bins.
The analysis of the south polar region (Schorghofer & Williams 2020) used six seasonal and 24 diurnal bins.However, this split between diurnal and seasonal bins turned out not to be appropriate for the north polar region, and the analysis here is based on two seasonal and 96 diurnal bins instead.Figure 1 shows the spatial distribution of the time coverage.The median of the time coverage is 167 out of 192 bins (87%).
In the south polar region, we find that the cold trap area determined from 96 × 2 versus 24 × 6 maps is very closely the same, but not so in the north polar region, which may be related to the distribution of the PSRs.In general, PSRs in the south are larger and closer to the pole, and the topographic relief is higher.Seasonality is less important in the north than in the south (Williams et al. 2019), so having fewer diurnal bins in favor of more seasonal bins may lead to unnecessary smoothing.Time-averaging requires uninterrupted time series, so missing data need to be interpolated.This interpolation is carried out in the frequency domain, without any spatial interpolation, using the same data processing pipeline we used for the south polar region, described in Schorghofer & Williams (2020).If data are missing in any of the 96 × 2 bins, they are filled with data from a neighboring λ s time bin.If that neighboring bin is also empty, a 4 × 2 block average of nearby subsolar longitudes and the two seasonal bins is used.Temperatures were slightly smoothed with a Fourier filter along the diurnal frequency axis, and the impact of this smoothing on the cold trap area is described below.
A map of mean temperatures for the north polar region, obtained by averaging the 192 maps, is shown in Figure 2.

Northern Cold Traps
Sublimation rates of ice in vacuum, E, are calculated based on the Hertz-Knudsen formula and the vapor pressure of ice, with quantitative details provided in Schorghofer & Williams (2020).A temperature of 110 K corresponds to a sublimation rate of 114 kg m −2 Ga −1 .The traditional criterion for a cold trap is that the maximum surface temperature, T max t ( ), is smaller than a threshold temperature T thres , or, equivalently, Here, a criterion for the time average of the sublimation rate, , is also used.The time average is the physically more relevant quantity.Figure 3 shows maps of the maximum potential sublimation rate, E max t ( ), and the time-averaged potential sublimation rate, mean t (E), in the north polar region.Figure 3(d) shows mean t (E) versus max t (E), which reveals that the time-averaged sublimation rate is always at least 10 times lower than the maximum sublimation.The ratio is limited by the number of time bins available (192), and it is at some locations undoubtedly even more extreme.
Figure 4 shows cold trap areas for various thresholds and methods.For the smoothed data sets, measured temperatures  instead of max t (E) results in larger areas of 7.1 × 10 3 km 2 , unsmoothed, and 7.7 × 10 3 km 2 , smoothed.The time average is the physically more accurate criterion, whereas the peak value is easier to evaluate.The cold trap area based on mean t (E) is 33% larger than for 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 Ga −1 .Figure 4 also shows the PSR area for the same geographic region and spatial resolution, which is significantly larger than the total cold trap area.
Individual cold traps were identified with a function that finds contiguous patches of pixels.Poleward of 80°N, there are 19 cold traps larger than 100 km 2 , 99 larger than 10 km 2 , and 327 larger than 1 km 2 , when based on a threshold of 100 kg m −2 Ga −1 for mean t (E) and the Fourier-filtered data set.

Comparison of PSRs and Cold Traps
PSRs are defined as areas that never receive direct sunlight, whereas cold traps are defined by a sublimation rate threshold.Areas permanently without direct sunlight might not be cold traps, because terrain irradiance contributes to the energy balance.The global PSR area is known to be far larger than the global cold trap area, even poleward of 80°latitude (Figure 4), so cold traps lie within PSRs.
However, this is not necessarily the case at each location.In principle, grazing sunlight can be of such low intensity that sublimation remains below the threshold (Rubanenko 2015).The direct solar flux Q that corresponds to a 110 K equilibrium temperature T eq is given by . With a bolometric albedo of A = 0.11 and an emissivity of ò = 0.99, the threshold flux is 9 W m −2 .This corresponds to a Sun elevation of 0.4°relative to the local slope.The incident rays could be even steeper when the Sun is partly blocked.Another hypothetical reason for sunlit cold traps is that direct illumination could be so brief that the surface does not have time to warm up to its equilibrium temperature.Below we will estimate how brief such an event would have to be.Moreover, the time average dilutes the sublimation rate of such a brief event.Grazing sunlight or the thermal inertia of the surface could, in principle, give rise to briefly sunlit cold traps.
Maps of PSRs are available based on ray-tracing calculations of LOLA (Lunar Orbiter Laser Altimeter) topographic data as described in Mazarico et al. (2011).Data product LPSR_65N_240M_201608 (LRO-L-LOLA-4-GDR-V1.0) is a map of the PSRs in the lunar north pole.The map resolution is 240 m px −1 , true at the pole in polar stereographic (spherical) projection.
Figure 5(a) shows contours of PSRs and contours of peak temperature in the north polar region.(If mean sublimation rates were used, the cold traps would be larger.)Many of the PSRs are not cold traps, particularly the smaller PSRs.Most cold traps lie within PSRs, as some PSRs can still receive sufficient thermal radiation or scattered light energy to raise their temperature above the volatile stability threshold.However, at a few locations the cold trap appears to extend beyond the PSR boundary, which is unexpected.These appear to be sunlit cold traps.The difference is particularly large at  Also computed is the maximum direct solar irradiance over the calendar year 2025 at 120 m px −1 , using the ray-tracing model by Mazarico et al. (2018).As mentioned above, incident fluxes above 9 Wm −2 result in equilibrium surface temperatures above 110 K from the direct flux alone.Figure 5(a) includes contours of a maximum direct solar flux of 10 Wm −2 , and these align closely with the PSR boundaries (0 Wm −2 ).In Lenard Crater this contour lies within the Diviner-based cold trap.Figure 5(c) reveals that direct irradiance reaches hundreds of Wm −2 in these supposedly perennially cold areas.Hence, the Diviner temperatures must be undersampled at this location, as direct illumination alone is sufficient to bring temperatures above the sublimation threshold.
Figure 6 shows a time series for a location that is a cold trap according to Diviner data, but sunlit according to illumination calculations.The Sun rises three times each Draconic year, but only a fraction of the solar disk emerges above the horizon.This location receives direct sunlight only 0.6% of the time, which explains why warm periods were missed by Diviner.The longest continuous period of direct illumination from 2000-2029 at this location is 29 hr.The angle between the incident Sun rays and the local slope reaches as high as 10 • , because the surface is sloped.The Sun rises at approximately east-northeast direction (far from noon), and the pole-facing slope is tilted toward that direction.Although only a fraction of the solar disk emerges above the horizon (Figure 6(d)), the incoming solar flux is substantial during these brief periods and radiative equilibrium temperatures reach over 200 K (Figure 6(c)).A slight warming is discernible in response to the illumination events (Figure 6(c)), which is possibly (but not necessarily) caused by these heating events.The time-averaged potential sublimation rate based on these temperatures is 0.02 kg m −2 a −1 , orders of magnitude too high to qualify as perennial cold trap.
It could be argued that Diviner bolometric temperatures might be biased cold when solar incidence angles are high (grazing illumination), because unresolved shadows are especially cold and large (Sefton-Nash et al. 2019).However, the brief periods of direct illumination were not sampled by Diviner, not only for the one year of data shown in Figure 6(c), but over the entire 15 Draconic years of Diviner observations, so, at least at this location, the discrepancy is definitively due to undersampling.
Terrain irradiance adds to the temperatures calculated based on direct illumination, whereas subsurface heat storage can delay the temperature increase.To estimate whether radiative equilibrium is reached during these brief periods of direct illumination, one can balance the incoming flux with the heat storage.Enough energy is available when Here Q abs is the absorbed solar flux, ρ density, c heat capacity, and Δz is the depth scale a heat pulse penetrates over the time Δt, which in turn can be estimated by z t k D » D , where κ is the thermal diffusivity.Combined, is the thermal inertia.By definition, the equilibrium temperature is given by T Q eq abs 4 s = . Hence, the approximate criterion becomes For T eq = 200 K and Δt ≈ 20 hr, the left-hand side is 120 J m −2 K −1 s −1/2 , whereas a typical nighttime thermal inertia for the lunar surface is 55 J m et al. 2017).Hence, these events are long and intense enough to raise surface temperatures close to the equilibrium temperature.
In the example shown in Figure 6, equilibrium temperatures are far above cold trap temperatures, so this location cannot be expected to act as perennial cold trap even when the inertia effect is taken into account.(It can be argued that if significant amounts of ice were present, the thermal inertia would be much higher and the peak temperature would be lower, but it is unclear whether this would suffice.) We conclude that the true cold trap area is smaller than can be determined from Diviner measurements, due to brief periods of direct sunlight that were not sampled by Diviner but inferred from ray-tracing calculations.Clipping areas found to be directly illuminated from the Diviner-based cold traps would be a step toward a more accurate map of cold traps, and might result in a significantly smaller cold trap area than determined here, especially close the north pole where Diviner coverage is low (Figure 1).

Subsurface Ice Stability
To evaluate potential subsurface ice stability, we modeled subsurface temperatures.The one-dimensional heat equation is solved with the time-varying Diviner temperatures as upper boundary condition for each of over 5 million pixels, as described in Schorghofer & Williams (2020).Depth-dependent thermal properties are assumed with the parameters given in Hayne et al. (2017) derived from the rock-free regolith temperatures of Bandfield et al. (2017).
One notion of "subsurface ice stability" is long-lived buried ice.Buried ice experiences a lower peak temperature than the surface.Moreover, the vapor has to overcome a diffusion barrier, so it experiences less sublimation loss than ice exposed on the surface.Sublimation rates are calculated based on the modeled subsurface temperatures and time-averaged over the last 12 diurnal cycles.These rates are then divided by a depthdependent factor that represents the diffusion barrier, as described in Schorghofer & Williams (2020).A threshold can then be applied to this three-dimensional map of sublimation rates to obtain a map of depth-below-the-surface for a specific sublimation threshold.
Figure 7 shows results for the sublimation rate of buried ice in the north polar region.Where subsurface ice is stable, the depth scale for the sublimation threshold is determined by the thermal skin depth.The shallow stability depths are explained by the small thermal skin depth.Any ice needs to be buried quickly to be subject to this reduced sublimation rate, so subsurface stability represents merely an upper bound for the occurrence of ice and only a small fraction of this area can be expected to actually harbor ice.
Another notion of subsurface ice stability is subsurface cold trapping.Water molecules that diffuse from the surface downward can be cold-trapped at depth.Downward diffusion is accelerated by diurnal temperature cycles ("vapor pumping effect").This process captures only a small fraction of water that was delivered to the polar regions, but over a large area (Schorghofer & Aharonson 2014;Reiss et al. 2021;Schorghofer 2022).
The same subsurface temperature calculations are used to evaluate subsurface cold trapping.If the time-averaged desorption rate on the surface exceeds the sublimation rate of potential ice at depth, subsurface cold trapping is expected to occur.The surface concentration of water molecules is determined by a differential equation that integrates the amount of adsorbed H 2 O over time, assuming a rate of influx and a rate of destruction by space weathering.Further details of this model are described in Schorghofer & Aharonson (2014).
Figure 8 shows the "pumping differential," which is the difference between the desorption rate on the surface and the sublimation rate in the subsurface assuming a gradual delivery of 1 nm a −1 of water to the polar region.Where this difference is positive, low concentrations of ice are expected to be coldtrapped in the shallow subsurface.The total area poleward of 80°N with a positive pumping potential is 109 × 10 3 km 2 , and the pumping differential is larger than half the supply rate over 94 × 10 3 km 2 .These areas are lower than was estimated in Schorghofer & Aharonson (2014) based on cruder temperature inputs, but 1 order of magnitude larger than the surface cold traps.These areas could be responsible for the excess hydrogen that has been detected outside of cold traps (e.g., Sanin et al. 2017;Gläser et al. 2021).

Conclusions
Diviner temperature measurements for the north polar region were binned according to subsolar longitude (local time) and subsolar latitude (summer/winter).Empty time bins were filled by frequency-domain interpolation.These temperature time series are then used to calculate time-averaged sublimation rates of potential surface and subsurface ice, and result in more accurate maps of potential ice reservoirs in the lunar polar regions.
The map of north polar cold traps, defined by long-term potential sublimation rate, is shown in Figure 3(b).The total cold trap area poleward of 80°N based on Diviner data at 240 m px −1 resolution is about 7 × 10 3 km 2 .Some Diviner-defined cold traps extend into sunlit areas, even when the peak sublimation rate is used for the threshold that defines cold traps (Figure 5), but models of direct illumination reveal that these areas are briefly illuminated each Draconic year, e.g., for one Earth day (Figure 6).Even after more than 15 Draconic years of near-continuous measurement, nadir-pointing Diviner data are not available at sufficiently high resolution in terms of subsolar and ecliptic longitudes.To reduce remaining uncertainties in the extent of cold traps, calibrated terrain irradiance models could be applied to existing high-resolution topographic maps.Subsurface temperatures were modeled using time-varying Diviner surface temperatures as upper boundary condition.Model-predicted maps for the loss rate of (potential) relic buried ice and (potential) subsurface cold trapping are shown in Figures 7 and 8.
Maps of cold traps, surface temperatures, potentially longlived buried ice, and subsurface cold traps are available online at Schorghofer & Williams (2024).

Figure 1 .
Figure 1.Time coverage of available data for 2009 July to 2024 January poleward of 80°N, defined as the percentage of the 96 × 2 = 192 temporal bins populated with Diviner data measurements.

Figure 2 .
Figure 2. Mean surface temperature poleward of 80°N averaged over diurnal and seasonal cycles.

Figure 3 .
Figure 3. (a) Maximum potential sublimation rate, E max t ( ), poleward of 80°N with logarithmic color scale.Black contours correspond to 100 kg m −2 Ga −1 . (b) Long-term average of the potential sublimation rate, mean t (E).(c) Comparison of contours for a 100 kg m −2 Ga −1 sublimation rate for E max t ( ) and mean t (E).The red contours are plotted on top of the blue contours, and always enclose a smaller area than the blue contours.The background map is a shaded topographic relief.(d) Correlation of E max t ( ) with mean t (E) on logarithmic scales.

Figure 4 .
Figure 4. Cold area poleward of latitude as a function of sublimation rate threshold and for various methods.The total PSR area based on a 240 m px −1 LOLA map is also shown.

Figure 5 .
Figure 5. (a) Comparison of PSRs and cold traps in a portion of the north polar region.Occasionally, Diviner-based cold traps extend beyond the PSR.The background map is a shaded topographic relief, and the north pole is at the lower right corner.(b) Fraction of diurnal-seasonal bins occupied with Diviner data.Contours are the same cold traps as in (a).(c) Maximum direct irradiance over the year 2025.

Figure 6 .
Figure 6.(a)-(c) Time series for the year 2019 at stereographic coordinates (−150 km, +62 km), the location indicated by a black star in panel (e) in Lenard Crater.This location briefly receives direct sunlight three times during the year, when the Sun partially rises above the horizon.(c) Radiative equilibrium temperatures are calculated based on the direct irradiance, assuming an albedo of 0.11 and an emissivity of 0.99, and compared with channel 9 Diviner measurements.(d) Elevation of the horizon and the Sun.Dots are spaced by 1 hr.The azimuth is measured clockwise relative to the right-pointing horizontal axis in map (e) and part of solar disk emerges above the horizon near summer solstice at sunrise.(e) Topography of the area relative to sphere.Contours are spaced 250 m, and the elevation varies between −3500 m and +1504 m within the region shown.The lowest horizon elevations, as viewed from the study location (black star), are on the far side of Hermite Crater (azimuths ∼30-80°), beyond the area shown.

Figure 7 .
Figure 7. Depth to loss rates of 1 kg m −2 Ga −1 poleward of 80°N.White indicates ice is stable at depths larger than 20 cm or not stable at any depth.