History of Ceres's Cold Traps Based on Refined Shape Models

Dwarf planet Ceres hosts permanently shadowed areas in its polar regions (Hayne & Aharonson 2015; Schorghofer et al. 2016; Ermakov et al. 2017; Platz et al. 2017), and these regions are an interesting analog to Mercury and the Moon (Lawrence 2017). Ceres's permanently shadowed regions (PSRs) were mapped using Framing Camera (FC) data from the Dawn spacecraft, and, thanks to scattered sunlight, bright deposits were discovered in a fraction of the PSRs (Schorghofer et al. 2016; Platz et al. 2017). The axis tilt (currently 402) can be robustly integrated backward in time and periodically oscillates between 2° and almost 20° with a period of 24.5 kyr (Bills & Scott 2017; Ermakov et al. 2017). The bright crater floor deposits (BCFDs) may correspond to the extent of PSRs in the past rather than at present day, which explains why most PSRs have no BCFDs. To arrive at a clearer understanding of the nature of cold-trapped ice deposits on Ceres, we construct improved shape models of PSR-hosting craters and revisit the relation between BCFDs and PSRs. We also model temperatures in the shadowed regions.

Before turning to stereophotoclinometry (SPC) models for individual craters, we consider results from the global stereophotogrammetry (SPG) model, which was not available in its current form at the time of previous ray-tracing studies (Schorghofer et al. 2016; Ermakov et al. 2017). The Ceres shape model DAWN-A-FC2-5-CERESHAMODTMSPG-V1.0 was constructed from High Altitude Mapping Orbit (HAMO) FC images through SPG (Roatsch et al. 2016) with a resolution of 136.7 m pixel⁻¹.

We conducted simulations to determine the locations in year-round shadow (PSR) over the range of possible Ceres obliquity values, from its minimum of ∼2° to its maximum of ∼20°. As the underlying shape model was obtained from present-day imagery with limited illumination around PSRs, it may not reliably reproduce the PSR area for larger-than-current obliquity values. The simulations were conducted over a ±330 km box around the north pole at the SPG model resolution, with corners at latitude ∼38° N and the side centers at latitude ∼50° N, so the spatial coverage is only complete poleward of 50° N. From each location, the surrounding topography's angular elevation was first computed from the shape model to determine the terrain horizon, following Mazarico et al. (2018). For each obliquity value, similarly to O'Brien & Byrne (2022), we compared these horizon masks to the maximum solar decl. over a Ceres year, assuming a circular orbit.

Figure 1 shows the areas that can be in year-round shadow, with color-coding indicating the maximum obliquity, above which the Sun will be visible at least in summer. The north polar PSR area derived from the current global SPG model in the current orbital configuration (4° obliquity) is 2282 km² when projected onto a sphere. This is greater than the 1800 km² obtained from ray-tracing with an early SPG model and a more restricted geographic domain (Schorghofer et al. 2016). Platz et al. (2017) determined the PSR area of the north polar region by stacking many images acquired near summer solstice and obtained 2130 km².

Zoom In Zoom Out Reset image size

The most equatorward PSRs are on the pole-facing interior slopes of craters. Within about 10° of the pole, PSRs tend to occupy crater floors. The largest PSRs are found in the craters Dikhan and Mlezi (Figure 1), but these, as with most other present-day PSRs, do not harbor any BCFDs. Table 1 lists the largest PSRs by area. In total, there are 63 PSRs with areas larger than 10 km², and these account for 77% of the total PSR area. (All areas are here projected on a sphere of 470 km.) Table 2 lists the northernmost (potentially coldest) and most equatorial PSRs. No PSR larger than 2 km² exists within 3° of the north pole. The most equatorward PSR of any size was identified at 49.1° N.

Figure 2 shows how the total area in permanent shadow in the north polar region varies with obliquity and with latitude. It illustrates how spatially limited PSRs are at high obliquity, as discussed in more detail in Section 4.

Zoom In Zoom Out Reset image size

SPC identifies control points in images by correlating illuminated topography and reconstructs short-scale digital topography and albedo maps from brightness-derived slope constraints. A small patch of topography, which we call a maplet, is constructed by solving for the slope and relative albedo at each grid point by minimizing the summed squared residuals between the image brightness for that maplet vertex and the predicted brightness over a large set of images. If a point is not illuminated in any of the images, its maplet vertex is left undefined. The slopes are then integrated, conditioned by external data such as heights from overlapping maplets, limbs, and differential stereo, to produce the final maplet.

Maplets here are 99 × 99 pixels large at 100 m pixel⁻¹ (about 10 ×10 km) or 40 m pixel⁻¹ (about 4 × 4 km). Ordinarily, SPC identifies shadows as areas with an image brightness below some preset threshold, and then these areas are masked. The Dawn cameras were sensitive enough to detect secondary illumination in the shadows of craters, so the threshold was adjusted, setting it to be 10% of the brightest part of the image extraction.

In order to use data from secondary illumination, an experimental procedure was developed whereby data from mainly the Low Altitude Mapping Orbit (LAMO) images were extracted from areas corresponding to the undefined area of the maplet, with a little spill allowed around the edges. These extracted data are stretched to produce a new image with the same spacecraft state as the original image but with an illumination sourced from an anti-Sun azimuth and an elevation of about 10°. Stretched images were added to the directly illuminated images to generate a new maplet. In standard SPC calculations, pixels that are in shadow are excluded (masked out). However, there is a wealth of information in the shadows of Ceres because of secondary illumination by light scattered by the sunlit portion of the crater. In this enhanced process, we use the same imagery but make adjustments so that SPC can use the dim pixels in the shadowed regions. First, we ingest a duplicate copy of each image. Then, the image brightness is stretched and terrain with primary illumination is masked out. Using these stretched images, we build a small digital elevation model over this region using only the data from the shadowed region. This maplet has the topography of only the shadowed region and no information about the directly illuminated regions.

Because SPC uses stereo and photoclinometry, a light source is needed to conduct the calculations. This is done by reversing the Sun vector and then setting the elevation to 10° above the horizon. The 10° angle corresponds to the approximate aspect ratio of a crater. We experimented with different Sun elevation angles, but that did not have a substantial effect on the outcome. The bright spots within the shadowed areas are too small to exhibit much internal structure and act primarily as stereo elements. HAMO data were included to improve the stereo solution. The shape solution depends on the relative angles, whereas the photoclinometric effect is comparatively small. For a similar reason, intrinsic brightness variations within the shadowed region (ice is brighter) will not change the outcome significantly.

Finally, we combine the maplet containing topography for illuminated regions with the maplet that contains the topography for the shadowed region. The small overlap aids in this process. Overall, this technique provides a consistent solution for the shape inside and outside the shadowed regions.

The Dawn spacecraft entered its mapping orbit around Ceres close to northern summer solstice in 2015. Although the spacecraft ultimately lasted until southern summer solstice over 2 yr later, it had ascended to a higher orbit by then, and the SPG and SPC shape models that can be constructed near the south pole do not have nearly the high resolution and quality available for the north pole. Hence, our analysis of the south polar region is limited to a single PSR, which was already identified in Ermakov et al. (2017).

4.1. Study Sites

4.2. PSRs Based on SPC Models for Individual Craters

The extent of bright deposits is compared with the model-derived extent of PSRs at various obliquities. Shadows were determined by ray-tracing, using python-flux (Potter et al. 2023), which in turn calls the Computational Geometry Algorithms Library (CGAL). Calculations were performed for circular orbits with 360 time steps over one solar day at summer solstice. Some of the FC images were manually aligned with the topography-derived maps, as the pointing information of the spacecraft is not always sufficiently accurate.

At Ceres's semimajor axis, the solar disk has a diameter of 11' and can thus be treated as a point source. The direct incoming solar irradiance is ${Q}_{\mathrm{direct}}=({S}_{\odot }/{R}^{2})\,v\cos \theta$, where S_☉ is the solar constant, R is the distance from the Sun, θ is the local incidence angle, and v is either 0 or 1, depending on whether the facet is shadowed or has a direct line of sight to the center of the solar disk.

Figure 3 shows mapped FC images, stretched in grayscale to focus on the shadowed portions, whereas the sunlit portions appear white. Superimposed on these images are color contours that represent the extent of PSRs for different values of the axis tilt.

Zoom In Zoom Out Reset image size

In Bilwis crater (NP1, #1), the BCFD is aligned with the PSR at 10° axis tilt. In Zatik crater (NP5, #2), the extent of the BCFD roughly follows that of the PSR at about 8° axis tilt. In crater NP7 (#3), the BCFD aligns with the current PSR, but the PSR contours are close to one another, allowing little differentiation between axis tilts. In unnamed crater #4, the bright deposit has an extent close to that of the PSR at 10° axis tilt. NP19 and NP26 were not considered to have BCFDs by Platz et al. (2017) but are included in Ermakov et al. (2017). In Damia crater (NP19), no PSR exists at 10°, but the slightly brighter area on the crater floor is very faint. In Enkimdu crater (NP26), the bright area on the crater floor is also faint and does not have a compelling BCFD, consistent with the lack of a PSR at 10° obliquity. The study site in the south polar region (SP1) has a bright deposit within a fainter bright deposit. The central BCFD best aligns with a PSR at about 4° axis tilt but allows little differentiation with axis tilt.

Unlike the previous results by Ermakov et al. (2017), these improved shape models suggest that no PSRs exist at an axis tilt of 20° and therefore these craters have no truly permanent PSRs. A PSR starts to emerge in Bilwis crater at about 18°, and they emerge at lower obliquities at the other six study sites. This implies that the ice deposits are remarkably young.

Figure 4 shows the recent history of the obliquity, as determined by backward integration in Ermakov et al. (2017). Currently, the obliquity increases from a recent minimum of 24 only 1.3 kyr ago. The most recent obliquity maximum of 185 occurred 13.9 kyr ago.

Zoom In Zoom Out Reset image size

The results in Figure 3 suggest that BCFDs formed at an axis tilt less than 10°. The BCFDs in NP7 and SP1 would be consistent with present-day PSRs but may be sensitive to remaining uncertainties in the shape models. An obliquity of 10° was reached 5.8 kyr ago (Figure 4). All BCFDs are less than 6 kyr old. An earlier delivery or exposure of ice would have resulted in larger or no BCFDs.

4.3. Can the Ice Deposits Endure Direct Sunlight?

To estimate the sublimation loss a bright deposit would experience when in direct sunlight, the one-dimensional heat equation is solved for a single flat and unobstructed location using the numerical model of Schörghofer (2023). The thermal inertia of the bright deposit could be as high as ∼2000 J m⁻² K⁻¹ s^−1/2 for consolidated ice, or much lower if in powder form. It is unlikely that ice on the surface would remain consolidated over a long time period, but calculations are carried out for a range of thermal inertias. The approximate latitude of Bilwis crater is chosen, the northernmost study site. The albedo was assumed to be twice the surface average of Ceres, based on the brightness ratios given in Ermakov et al. (2017). Time integration proceeded in 50 steps per solar day, for each of 4446 sols along one orbit. The thermal model was equilibrated over 20 orbits, and peak and average temperatures and sublimation rates were calculated for the last orbit.

Figure 5 shows the results for an obliquity of 18°, the approximate value for the most recent obliquity maximum. As long as the thermal inertia is below about 1000 J m⁻² K⁻¹ s^−1/2, potential sublimation rates are above 10 mm yr⁻¹. Over the 6 kyr before all of the bright deposit became continuously shadowed, meters of ice would have been lost while exposed to direct sunlight. Such massive losses could only be sustained for a short time. Higher albedo would reduce the sublimation loss (Hayne & Aharonson 2015), but very high albedo values cannot be sustained in the long term owing to darkening by space weathering, nor are they observed. The observed bright deposits disappear when exposed to direct sunlight at high obliquity and must therefore be young.

Zoom In Zoom Out Reset image size

Whereas PSRs are defined by shadows, cold traps are defined by temperature, or, more precisely, by the (potential) sublimation of ice into vacuum, which is determined by temperature. That threshold is conventionally set at a peak temperature of about 110 K, which corresponds to a sublimation rate of 0.15 m Gyr⁻¹ of water ice. At 120 K, the sublimation rate is already two orders of magnitude higher, 14 m Gyr⁻¹, which for Ceres is still only 0.3 mm per obliquity cycle.

5.1. Temperature Model

Temperature modeling is also carried out with python-flux (Potter et al. 2023), the same software that was used for the shadow modeling. A view factor matrix F is constructed (e.g., Lagerros 1997), which is then used to scatter radiation between surface elements (terrain irradiance). The dimensions of the shape models (19,208 facets each) are small enough to allow direct computation of these view factor matrices. Since the albedo is small, most of the terrain irradiance is in the infrared rather than the visible spectrum.

Assumed parameters are S_☉ = 1365 W m⁻² for the solar constant, a Bond albedo of A = 0.04 (Li et al. 2019), and an infrared emissivity of ε = 0.95. Two radiosity equations are solved, one for visible wavelengths, Q_vis, and another for infrared wavelengths, Q_IR,

$$\begin{eqnarray}&&(1-{\boldsymbol{F}}A){Q}_{\mathrm{vis}}={Q}_{\mathrm{direct}}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&(1-{\boldsymbol{F}}){Q}_{\mathrm{IR}}={\boldsymbol{F}}(1-A){Q}_{\mathrm{vis}}.\end{eqnarray} \tag{ 2 }$$

The total absorbed irradiance is

$$\begin{eqnarray}&&(1-A){Q}_{\mathrm{vis}}+\varepsilon {Q}_{\mathrm{IR}}=\varepsilon \sigma {T}_{\mathrm{eq}}^{4},\end{eqnarray} \tag{ 3 }$$

where σ is the Stefan–Boltzmann constant. This provides the equilibrium temperature T_eq.

A circular orbit was assumed. Equilibrium temperatures were calculated for 180 time steps over the solar day at summer solstice to determine the peak temperature. The thermal inertia of Ceres is small, 40–160 J m⁻² K⁻¹ s^−1/2 according to measurements at millimeter wavelengths (Li et al. 2020), so surface temperatures should quickly approach the equilibrium temperature. Cohesive ice increases the thermal inertia, so the model temperatures should be considered to be representative of ice-free soil.

5.2. Cold Traps for Supervolatiles

Figure 6 shows contours of the maximum equilibrium temperature at each of the seven study sites for an axis tilt of 4°. The area that remains below 80 K over the entire orbit is negligibly small at current axis tilt (<1 km² for all study sites combined). Cerean PSRs are surprisingly warm considering that some lunar PSRs have temperatures well below 80 K (Williams et al. 2019). This is due to the larger axis tilt (4° instead of 15) and the absence of PSRs very close to the north pole of Ceres (Figure 1). The study site closest to the pole is Bilwis crater (862 N), where the Sun rises almost 8° high at summer solstice, leading to considerable indirect illumination. These surface temperatures are consistent with results for idealized bowl-shaped craters on Ceres at the same Sun elevation (Schorghofer et al. 2016).

Zoom In Zoom Out Reset image size

The model temperatures provide insight into whether any supervolatiles (also called hypervolatiles) could be cold-trapped on Ceres. Zhang & Paige (2009) list threshold temperatures for numerous species, including 107 K for H₂O, 80 K for NH₃ · H₂O (ammonia clathrate), 66 K for NH₃, 62 K for SO₂, and 54 K for CO₂. The precise temperature threshold for each species depends on the chosen sublimation loss threshold. Barely any cold trap area is available even for NH₃ · H₂O.

Bilwis crater (NP4, #1) hosts the northernmost large PSR (Table 2), so it is unlikely that any significant area capable of trapping supervolatiles exists in the north polar region. Supervolatiles can only be trapped seasonally. Temperature calculations for an axis tilt of 2° (not shown), when the PSRs are coldest, also reveal little area that remains below 80 K over the entire orbit (<3 km² for all study sites combined, with most in NP19).

For the current eccentric orbit, the incoming flux is 16% lower at perihelion than for a circular orbit, so equilibrium temperatures could be lower by up to one-quarter of this percentage, 4% of the absolute temperature. This is not enough to change the conclusion about the absence of supervolatile cold traps. In any case, the orbit and spin axis precess over time.

5.3. Sublimation of Water Ice

The sublimation rate of ice into vacuum is a function of temperature and can be calculated based on the saturation vapor pressure of ice and the Hertz–Knudsen relation. Equilibrium temperatures were calculated every 15° along the orbit (including the day of summer solstice) for one solar day each, using 90 time steps for each sol. Sublimation rates for these temperatures are then time-averaged over each sol and over one orbit.

Figure 7 shows sublimation rates for Bilwis crater. In this crater, PSR contours are spread out spatially (Figure 3), so it provides a better constraint than other study sites. The BCFD best aligns with the PSR at an axis tilt of 10°. Sublimation rates over the bright deposit for an axis tilt of 10° (5.8 kyr ago) are <10⁻⁴ mm yr⁻¹. Since then, the axis tilt decreased to its current ∼4°, and sublimation rates were thus even lower. Over the last 5.8 kyr when the BCFD was continually shadowed, the total sublimation loss is on the order of 0.1 mm. For the other study sites, PSR contours also remain inside the 10⁻⁴ mm yr⁻¹ contour, so similar integrated sublimation rates apply.

Zoom In Zoom Out Reset image size

Outside of PSRs, sublimation rates quickly exceed 1 mm yr⁻¹ (Figure 7). Hence, the faint bright region that surrounds the BCFD is unlikely to be water ice, as any ice would quickly be lost even in the short time since the last obliquity minimum. The faint bright region largely corresponds to the flat crater floor.

5.4. Sunlit Ice in Zatik Crater

At Zatik crater, the bright deposit extends beyond the PSR. The sunlit portion allowed Platz et al. (2017) to spectroscopically identify the bright deposit as water ice.

Figure 8 shows an FC image of Zatik crater stretched so that the sunlit ice west of the shadow is visible, whereas the majority of the BCFD lies hidden in the dark shadow. Topographic contours reveal that the sunlit ice deposit is located at a steep portion of the crater. Slopes go above 45° in this area, well above a typical angle of repose. The sublimation rates reveal that the sunlit portion is not only outside the PSR but also outside the cold trap area, with sublimation rates above 1 m yr⁻¹. Due to the high albedo and thermal inertia of the ice, these model sublimation rates may be overestimated, but sublimation is undoubtedly rapid. This portion of ice is outside even the maximum extent of PSRs at 2° axis tilt (Figure 3) and could not have sustained high loss rates for thousands of years. Only a recent exposure is consistent with the presence of ice at this sunlit location. A very recent mass wasting event may have exposed ice, as has been observed outside the polar regions (Combe et al. 2016, 2019). Ermakov et al. (2017) pointed out that the deposit in Zatik crater is brighter than the other BCFDs, which also points toward a young age.

Zoom In Zoom Out Reset image size

The bright ice deposits in the PSRs could have several possible origins:

• 1.  
  Ice could be cold-trapped from the water exosphere. The gravity of Ceres is strong enough that some water molecules that travel at thermal speed can return to the surface. The global ice mantle that gradually retreats by sublimation and diffusion can supply such an exosphere (Tu et al. 2014; Schorghofer et al. 2016; Landis et al. 2017; Schorghofer et al. 2017; Landis et al. 2019).
• 2.  
  An impact into the global ice-rich crust of Ceres (Prettyman et al. 2017) could have produced a short-lived atmosphere, and the cold traps would have then captured a fraction of this water atmosphere.
• 3.  
  Ice could have been exposed by mass wasting. In the polar regions the depth to ice is likely only centimeters or less (Prettyman et al. 2017), so even small amounts of erosion could expose ice on the surface. The few ice deposits that have been detected spectroscopically outside the polar regions are indeed often associated with landslides (Combe et al. 2016, 2019), and the sunlit portion of the ice deposit in Zatik crater is best explained by a recent mass wasting event. Ice would need to have been exposed or deposited recently over all areas where BCFDs reside. Statistically, events this recent must have occurred frequently, and over time this process would be expected to result in significant changes to the crater topography.
• 4.  
  Ice could be exposed by many small impacts. Any exposed ice would rapidly sublimate except in the cold traps. The uniform appearance of the BCFDs, their young age, and the synchronous appearance in multiple craters make this explanation unlikely.

The possibility best aligned with the observations is cold trapping from a transient water atmosphere, although mass wasting most likely played a role at Zatik crater and may be a viable explanation for the origin of other BCFDs as well. The polar ice deposits on Mercury are also most readily explained by an impact event that formed a transient water atmosphere and subsequent cold trapping (Chabot et al. 2018; Ernst et al. 2018).

Another hypothesis would be that the BCFDs do not consist of water ice at all. However, the strong relation with shadows and the spectroscopic identification of one BCFD as water ice are difficult to challenge. Moreover, Ceres has an abundance of water ice. The ice deposits in the Cerean PSRs indicate an active water cycle; ice is either repeatedly captured and lost or frequently exposed, or both.

The cumulative area of PSRs in the north polar region of Ceres is about 2200 km² at present, based on the global SPG model. A PSR as far equatorward as 49° N has been identified. No PSR larger than 2 km² exists within 3° of the north pole.

We devised a technique to construct improved digital elevation models that uses secondary illumination (i.e., light reflected from the crater walls). A separate topographic model for the shadowed region is constructed using a point source from the anti-Sun direction and then merged with the topography for the directly illuminated region obtained using standard SPC.

Based on seven study sites, no PSRs survive at the maximum axis tilt of 20°, and the BCFDs (the ice deposits) must therefore be recent (≲14 ka and likely ≲6 ka). If PSRs still existed at maximum axis tilt, the BCFDs could have accumulated over billions of years instead. The remarkably young age of the ice deposits exposed in the shadowed craters implies that delivery or exposure of ice is frequent.

Most BCFDs may align with PSRs at about 10° axis tilt (last reached about 6000 yr ago), which suggests that the delivery was due to a short-lived event rather than continuous. A recent small impact into the ice-rich crust of Ceres could have produced a short-lived atmosphere, and the cold traps could have then captured a fraction of this water atmosphere. Alternatively, avalanches have exposed ice and covered the crater floor with ice-rich material.

A terrain irradiance model was used to calculate equilibrium surface temperatures inside permanently shadowed craters. The lowest peak temperatures within PSRs are about 80 K, too warm for any common supervolatiles to be stored. These unexpectedly high temperatures are due to the lack of PSRs very close to the pole and the higher axis tilt compared to Earth's Moon.

We thank Anton Ermakov, Eric Palmer, and Jennifer Scully for insightful discussions. This material is based on work supported by the National Aeronautics and Space Administration under grant No. 80NSSC21K1033 issued through the Discovery Data Analysis Program (DDAP).

Maps of PSRs and maximum equilibrium temperatures of the study sites are publicly available at Schorghofer et al. (2023). The digital elevation models will also be available on the PDS Small Bodies Node.

Software: Computational Geometry Algorithms Library (CGAL, https://www.cgal.org), GMT (https://www.generic-mapping-tools.org), Octave (https://octave.org), Python (https://www.python.org), python-flux (https://github.com/sampotter/python-flux), Small Bodies Image browser (https://sbn.psi.edu/pds/sbib/).