Using Boulder Tracks as a Tool to Understand the Bearing Capacity of Permanently Shadowed Regions of the Moon

Permanently shadowed regions (PSRs) are abundant at the lunar poles. They experience no direct sunlight and reach temperatures as low as 30 K. PSRs are of interest as evidence suggests that some may contain water ice (H2O/OH‐), which could provide a record of the evolution of volatiles in the inner solar system. This water ice is also a critical resource for life‐support systems and rocket propellant. A better understanding of mechanical properties of PSR regolith, such as its bearing capacity, will help optimize the design of future exploration rovers and landers. Thirteen boulder tracks were identified on the edge of, or inside, south polar lunar PSR enhanced imagery and used to estimate the strength of the PSR regolith at latitudes of 70° to 76° in sites with maximum annual temperatures of 65 to 210 K. PSR boulder track features are similar to those observed in highland, mare, and pyroclastic regions of the Moon, implying similar properties of the regolith. Measured features were used to estimate bearing capacity for PSR regolith at depths of ~0.28 to 4.68 m. Estimated bearing capacity values suggest that these PSRs may be somewhat stronger than highland and mare regions at depths of 0.28 to 1.00 m. Bearing capacity in these PSRs is statistically the same as those in other regions of the Moon at depths of 1.00 to 2.00 m. The results of this study can be used to infer bearing capacity as one measure for the trafficability of lower‐latitude PSRs of the type measured here.


Introduction
The International Space Exploration Coordination Group (ISECG) Global Exploration Roadmap (2018) outlines the shared goals of the ISECG space agencies to expand human presence in low Earth orbit, to the Moon, and finally to Mars. A critical step in achieving these goals is continued exploration of the Moon and identifying and harvesting its available resources. Of particular interest in this endeavor are the permanently shadowed regions (PSRs) at the lunar poles. These regions are considerably colder than the rest of the lunar surface, reaching temperatures below 30 K . As a result, it has been theorized that these PSRs could trap and accumulate stable water ice, as well as other volatiles, and retain them for hundreds of millions of years (Siegler et al., 2018;Watson et al., 1961). Indirect evidence and more recent direct evidence support these findings using data from the Lunar Crater Observation and Sensing Satellite (LCROSS) lunar impactor and the Moon Mineralogy Mapper, respectively Li et al., 2018). These volatiles have a range of important applications for in situ resource utilization (ISRU), such as life-support systems and rocket propellant (e.g., Crawford, 2015). PSRs are also a valuable scientific target. The National Research Council outlined and ranked a series of scientific goals that should be considered when planning missions to the Moon (National Research Council, 2007). A future mission to a PSR would primarily address goals 4a, 4d, 7b, 7c, and 7d. These objectives pertain to the presence and behavior of volatiles on the Moon and the qualities of the regolith in which they reside. A robotic or human mission taking in situ measurements will be required to confirm volatile concentrations in the PSRs and assess the utility of any deposits for ISRU.
The success of a robotic mission is dependent on the ability of a rover to traverse the terrain. Thus, before planning a mission to a PSR, it is necessary to better understand the mechanical behavior of the regolith within these areas. Issues regarding extreme temperature lows, lack of solar energy, and static discharge (Jackson et al., 2011) in these regions are just a few of many concerns. Of particular concern for engineers is the soil behavior of the regolith in PSRs. The porosity of the regolith would influence the strength of the terrain and, thus, affect rover mobility. The porosity of PSR regolith has been inferred from the results of the LCROSS impactor where the resulting ejecta angles and flashes indicated highly porous material of 70% to equivalent depths of 2 to 3 m . Meanwhile, far-ultraviolet (FUV) reflectance properties of lunar PSRs showed FUV darkness, which could be a consequence of higher porosity,~70%, in the surface soil compared to regions outside PSRs with porosities of~40% . LRO Diviner indicates low thermal inertia values in the lunar high-latitude regions, which suggest that the upper few centimeters of regolith are likely highly porous (Hayne et al., 2017). Laboratory studies have also shown that low thermal cycling in PSRs could result in higher-porosity regolith (Metzger et al., 2018). As a result of these studies, one could assume that PSRs would be more difficult to traverse. If rovers and landers experience more sinkage into the regolith, becoming stuck could mean the end of a mission.
As it stands, the physical properties of the regolith in lunar polar PSRs are relatively unknown as there have been no in situ measurements taken to date. Experiments on lunar soil simulants have provided estimates of how the presence of water ice increases the soil strength (Gertsch et al., 2006;Pitcher et al., 2016). Gertsch et al. (2006) found that increasing water ice content in the simulant resulted in a stronger material, behaving like weak shale or mudstone at only 0.6% to 1.5% ice content. Similarly, the strength of icy rock has been studied and shown to be stronger than dry rock at equivalent temperatures (Atkinson et al., 2018). However, this experimental work did not provide estimates of trafficability. The following work will represent the first comprehensive study of trafficability in PSR regolith using data derived directly from lunar surface features.

Methodology
This study investigates the geomechanical properties of regolith within and on the edge of south polar lunar PSRs using rockfalls and their associated boulder tracks, both of which are abundant on the Moon Kumar et al., 2016). A south polar PSR map generated by McGovern et al. (2013) was used to select images, taken by the Narrow Angle Camera (NAC) on board the Lunar Reconnaissance Orbiter (LRO), of PSRs that may host a rockfall. South polar images were processed to enhance the visibility of features within the PSRs so that boulders and their tracks could be measured. Those measured values were then used to estimate bearing capacities using the methodology of Bickel et al. (2019).

Image Processing and Digital Elevation Model Processing
It is not immediately evident that boulders and boulder tracks can be discerned within shadowed regions. However, we discovered that features within PSRs can be observed using secondary light diffusely scattered from an illuminated portion of a crater containing a PSR. Adjusting the brightness and contrast of Lunar Reconnaissance Orbiter Camera (LROC) NAC images, followed by image filtering, illuminates those features sufficiently well to be measured. Image filtering is necessary because raw NAC images are affected by semisystematic vertical and horizontal stripes as well as white noise. The noise was exasperated by contrast enhancement and greatly impedes accurate measurements of boulder tracks within these areas. To suppress or remove the noise content, two separate filtering methods were implemented. We outline those methods here, while referring readers to detailed workflows and filtering information in Figure S1 and Tables S1 and S2 in the supporting information. A total of 33 LROC NAC images were downloaded and processed by either of the following methods.
The first image processing method forced an exaggeration of the dominant vertical noise in echocorrected NAC images using a sliding rectangular kernel. This technique was designed to detect and mask vertical lines with a width of a single pixel. Subsequently, a Standard Hough Transform (SHT; Hough, 1962) was applied to track the lines in the image and denote their beginning and end by using the SHT matrix. Next, these marked lines were removed from the image, and the resulting gaps in the image were interpolated using the values of the neighboring pixels. To minimize interference between adjacent lines, the filter was only run for lines that were not adjacent to another line. Once the lines were removed and interpolated, the filter was run again. This process was repeated for a specified number of times that was dependent on the noise level in the image. The resulting images were further refined with a two-pass 2-D median filter using a 3 × 3 pixel kernel to suppress the remaining white noise content. The filtered images were up-sampled using a bicubic interpolation in a 4 × 4 neighborhood to further increase image quality.
The second method was performed in ISIS3.5.1 and applied a low-pass filter kernel (1× vertical image dimension) to the echo-corrected NAC images to calculate a global column average of the images. The resulting low-frequency image was then high-pass filtered (1× horizontal image dimension), which produced a global albedo-suppressed average noise image. This image was then subtracted from the original input. Subsequently, a standard 2-D median filter within a sliding 3 × 3 pixel kernel was used to suppress the remaining white noise content.
Although all applied filter routines were nonaggressive, a small amount of local artifacts were observed (en échelon and cone-shaped noise patches), but they did not interfere with the performed measurements. Of these two methods, both filters were used on a case-by-case basis depending on which method produced a superior view into the PSR. Lastly, all optimized images were warped to a planetocentric polar stereographic projection and translated and exported for utilization in ArcGIS10.1. An example of a raw NAC image of a PSR and the resultant filtered image is shown in Figure 1.
Depending on local illumination conditions and spatial resolution of a NAC raw image, the quality of the filtered products was sufficient to produce digital elevation models (DEMs) within PSRs with spatial resolutions of 2 to 3 m/pixel using the Ames Stereo Pipeline (Shean et al., 2016). Details about the produced DEMs can be found in Table S3.

Identifying Boulder Tracks
In order to see boulder tracks within a PSR, there must be some secondary illumination of the PSR when the NAC images were obtained. At latitudes~70°to 76°, there is sufficient secondary illumination of craters so that processing of the high-resolution LRO NAC images (~0.5 m/pixel) reveals boulders and their associated tracks within the PSRs ( Figure 2). Although higher-latitude PSRs have been imaged (Haruyama et al., 2008) using data from Selene's Terrain Camera and Multispectral Imager, both those instruments have spatial resolution well below that of LRO's NAC. For example, the Terrain Camera imagery used to view the floor of Shackleton crater had a resolution of 10 by 10 m, far too coarse to measure boulder shape, boulder track width, and boulder track depth. A complete survey of that geographic band was performed so that all measurable boulder tracks were identified. All boulders identified as "in shadow" in NAC images were then compared to the PSR map produced by McGovern et al. (2013), derived from ray tracing techniques applied to a topographic data set from the LRO, to confirm that the boulders lie within a PSR or on a boundary. Ultimately 11 boulder tracks were measured inside PSRs, with two boulder tracks being identified as seasonally shadowed regions, SSRs (Kloos et al., 2019). In this work, SSRs are defined as areas of shadow, which are extensions of PSR shadows, but do not fall within the PSR boundaries identified by McGovern et al. (2013). An example of a boulder track identified inside a PSR is shown in Figure 3. We recognize that secondary illumination has been mapped in PSRs closer to the pole (Cisneros et al., 2017;Koeber & Robinson, 2013;Koeber et al., 2014); however, the spatial resolution of the images has to be artificially reduced to account for the blurring introduced by the satellite's motion in combination with longer camera exposure times. Consequently, this method would not be suitable for higher-latitude PSRs with current orbital imagers.

Calculating Bearing Capacity
Measurements of the boulders and their associated tracks are used to estimate ultimate bearing capacity. Bickel et al. (2019) compared two methods for estimating ultimate bearing capacity based on work by  (per McGovern et al., 2013) and the location of identified PSR boulder tracks in green. The blue band area represents the region where boulder tracks were observed in PSRs, located between~70°and 76°latitude, that is, the region where secondary illumination and image processing still allow to use NAC images without compromising spatial resolution. Terzaghi (1951) and Hansen (1970). An analytical solution provided by Hansen (1970) for shallow foundations provides a better and more realistic estimate of bearing capacity, because it allows for the adaptation of the bearing capacity estimation to the local slope and boulder shape . Here we apply the Hansen (1970) equation for bearing capacity, q f , with the same assumptions used by Bickel et al. (2019) for other regions of the Moon: where c is the cohesion of the soil; q 0 is the vertical stress within the soil; γ is the unit weight of the soil; B is the width of footing; N (c,q,γ) are the bearing capacity factors, which are derived from the angle of internal friction, φ; and d (c,q,γ) are the depth factors, s (c,q,γ) are the shape factors, g (c,q,γ) are the local slope inclination factors, and i (c,q,γ) and b (c,q,γ) are the load and foundation inclination factors, respectively. The unit weight of the soil is defined as ρ·g, where ρ is the density of the soil and g is the gravity of the Moon. For a detailed explanation of the Hansen (1970) equation and its assumptions, see Bickel et al. (2019).
Density, cohesion, and internal friction angle are essential properties for understanding the behavior of the lunar regolith. These properties have been estimated through several studies; however, the values are highly varied in literature (Carrier et al., 1991;. Calculated bearing capacities in other works vary depending on the values pulled from the literature (Hovland & Mitchell, 1969;Mitchell et al., 1973;Moore, 1970). Using the same conservative assumptions as Bickel et al. (2019), a cohesion value of 1 kPa  was used for this work. As in situ measurements have yet to be taken for PSRs, highland bulk density values of 1.39 and 1.66 g/cm 3 were applied for shallow (<30 cm) and deep (>30 cm) PSR soil, respectively , as PSRs occur in highland terranes with feldspathic composition and relatively low iron abundances (Lawrence et al., 2002;Spudis et al., 2013). The angle of internal friction (29°) was estimated from granular flow features at their static angle of repose on slopes within PSRs evaluated with high-resolution DEMs. For details on the production of those DEMs from filtered NAC images, see Table S3; for details about the determination of the static angle of repose, see Bickel et al. (2019).

Measurements
From the 33 NAC images that were processed, a total of 13 boulder tracks were identified in or on the edge of 5 PSRs that were found to be suitable for measurements needed to calculate the ultimate bearing capacity q f (see Table S4). Measurements were made on the long L and short B diameters of the boulder, the width of the associated track W, the length of shadows cast within the track S, and the local slope α. All measurements were made near the boulder responsible for the track (subscript b) and near the beginning of the track (subscript s). Taking measurements at the start and end of the track allowed for the approximation of a static system, as the boulders will have little kinetic energy . Thus, in general, each track produced two data points for analysis. The width and shadow length measurements were repeated three times for each of the individual tracks whenever possible in order to calculate an average and to minimize the influence of variations in track appearance, a consequence of irregular boulder shapes forming tracks of changing widths. Shadow length and track width measurements were performed in the same area of the boulder track. Accurate measurements of the shadow length required determination of the illumination direction; the measurement was then made parallel to that direction. The measurement was most accurate when the shadow length was perpendicular to the track .
Generally, the depth of a boulder track D was derived from the shadow length L and the incidence angle i of the Sun where Determining the track depth in a PSR is challenging as the incidence angle of sunlight is not directly responsible for the length of shadows; instead, the shadow is produced from light reflected from an illuminated face of a crater. If specular reflection was wholly occurring (Figure 4a), sunlight should be reflected out of the crater. As some sunlight is clearly being reflected into the PSR, diffuse reflection is occurring (Figure 4b), likely as a result of surface roughness on these slopes. Assuming that diffuse reflection results in a broad beam of light from the illuminated crater wall into the PSR, the effective incidence angle θ can be derived (Figure 4c). The illumination source was taken as the midpoint of the illuminated section of the crater wall H/2, which was derived from a 100-m resolution LRO Wide Angle Camera (WAC) stereo-derived DEM. Using the elevation difference between this reflecting surface and the location of the boulder h, as well as the horizontal distance between these points x, the effective incidence angle can be calculated: The boulder track depth D in PSRs can therefore be calculated as follows:

Results
We began with a qualitative study of PSR boulder tracks in NAC images. Boulder tracks that crossed the boundary from sunlit regions into PSRs (e.g., M1117841678LE/RE) show no significant difference in track morphology, suggesting similar geomechanical properties ( Figure 5). Additionally, when compared to the tracks identified by Bickel et al. (2019) in highland, mare, and pyroclastic regions (Figure 6), the tracks do not show any significant difference in morphology. Boulders rolling into PSRs penetrate sufficiently deep to leave measurable tracks but do not completely sink into the regolith as might be expected if the regolith was unusually porous and weak. If there is a large density contrast between sunlit and PSRs, then the track size would be significantly different. Yet no difference is observed.

Journal of Geophysical Research: Planets
Once all measurements were obtained, Hansen's bearing capacity formula (equation (1)) was applied. Estimated bearing capacity values for each terrane (PSRs, highland, mare, and pyroclastics) were plotted with respect to the slope angle where the measurements were taken. The data are presented as distinct ranges for depth of track (≤1 m, 1 to 2 m, and 2 to 6 m) comparable to those of Bickel et al. (2019). A least squares fit was applied with a shaded error bar derived from an estimate of the standard deviation of the error (Figure 7). Potential slope limits are also shown. Bearing capacity tends to decrease as the slope angle increases, a result of the reduced soil volume bearing the boulder (Castelli & Lentini, 2012;Meyerhof, 1957) and potentially a reduction of physical soil strength due to slow downslope movement Carrier et al., 1991). In addition, the bearing capacity is shown to be higher in deeper material, which would be expected as the vertical stress and bearing soil volume increases over depth. Generally all the regions plotted follow a similar trend.
A plot of bearing capacity with depth is shown in Figure 8 for PSRs compared to highland, mare, and pyroclastic data from Bickel et al. (2019). As in Figure 7, Figure 8 also shows that the bearing capacity of the regolith appears to increase with depth, which is consistent with Apollo data . Calculated bearing capacities increase with depth  for all terranes, with PSR and pyroclastic regolith being generally somewhat stronger than highland and mare regolith at equivalent depths. However, there also appears to be a lot of scatter in the bearing capacities and overlap of values among terranes, which might be connected to varying slope angles at the respective measurement locations.
The trends produced in Figure 7 were used to estimate bearing capacity for each associated depth range as a function of slope (Table 1). The uppermost regolith that has been measured with this technique in PSRs (0.28 to 1.00 m) appears to have higher average bearing capacities than highland and mare material at ≤1 m but ranges overlap. The greatest variability in bearing capacity at shallow depths lies with PSRs, while at greater depths the greatest variability lies with pyroclastics. The data obtained from the boulder tracks in SSRs fit the same trend as those obtained from PSRs, with no noticeable anomalies, indicating that the soil strength properties in the SSRs and PSRs measured at these 70-76°latitude locations are similar.
A statistical analysis of the spread of calculated bearing capacities over depth was performed to determine if there are significant differences between terranes (see Table S5). A two-sample T test assuming unequal variances (two-tail) was applied to the calculated PSR values from this work and each data set from highland, mare, and pyroclastic calculations from Bickel et al. (2019) for the ≤1and 1 to 2-m depth range providing p values based on a 95% confidence interval. The 2-to 6-m depth range could not be considered as there were too few data points for PSRs at this depth. Two analyses were performed, the first considered all bearing capacity values in the depth range of ≤1 and 1 to 2 m across all slope angles. However, as there are no data for relatively shallow slope angles in the ≤1-m depth range in the highland and mare regions, a second T test was performed for bearing capacity data that are obtained at slope angles of >15°only for the ≤1-m depth range. This allowed for a comparison of shallow data without a bias associated with slope angle dependence.
The analysis showed that for all data in the depth range of ≤1 m, the estimated bearing capacity in PSRs, as determined by the distribution of data, is statistically similar to lunar pyroclastic deposits, and statistically stronger than both highland and mare regions. However, the analysis performed in the >15°slope angle range indicated that PSRs are statistically stronger than pyroclastic deposits, highland, and mare regolith. Meanwhile, at the depth range of 1 to 2 m, PSRs have statistically similar bearing capacities to pyroclastic deposits, highland, and mare regions when considering all slope angles and slope angles of >15°only. In contrast, there is some indication that the bearing capacity is weaker at greater depths. For example, an estimated value of 271 kN/m 2 at 2 to 6 m in a PSR is lower than 314 kN/m 2 in highlands and 336 kN/m 2 in pyroclastic material at that same depth. We note, however, that there are not sufficient measurements from those depths to be statistically significant. It is important to note that there is considerable variation in any single terrane at all depths.

Likelihood of Ice in Lower-Latitude PSRs
Models indicate that water ice will be stable and form on the surface of PSRs at temperatures below~100 to 112 K (Hayne et al., 2015;Paige et al., 1992;Siegler et al., 2015;Zhang & Paige, 2009. To understand whether water ice could be present in the PSRs studied here, a maximum (~240 m/pixel) temperature map derived from the Diviner Radiometer  Note. The highland, mare, and pyroclastic data were obtained from Bickel et al. (2019). Data were omitted where the estimated bearing capacity values fall below 0 kN/m 2 or where a trend line with uncertainty could not be obtained (≤3 data points). Due to Hansen's (1970) assumptions, estimated bearing capacity approaches zero when local slope angles approach the angle of internal friction; thus, data are sparse in the 30°b in.

10.1029/2019JE006157
Journal of Geophysical Research: Planets data  was evaluated. The PSRs in this study have maximum surface temperatures of 65 to 210 K (see Table S6). Four of the boulders studied in this work, distributed between two PSRs, are found in locations that have temperatures that would be suitable for surficial surface ice to be stable, although it is not observed. The presence of CO 2 ice is unlikely in the PSRs studied here as they require temperatures of <55 K to be stable (Hayne et al., 2019;Zhang & Paige, 2009) and, therefore, these "warmer" PSRs may be more desirable targets for those wishing to obtain water ice without CO 2 impurities. A lunar polar ice map derived from Moon Mineralogy Mapper reflectance spectra (Li et al., 2018) suggests that surficial water ice (within a few millimeters of the surface) is not present at any of the studied boulder locations, suggesting that if ice has formed in the past, it would likely have redistributed deeper as a result of gardening of the upper regolith (Hurley et al., 2012). It should be noted that the PSRs studied here do not represent the PSRs found at higher latitudes that do have evidence to suggest surficial water ice.
The absence of surficial ice does not preclude the presence of subsurface ice. Water can thermally migrate through the regolith where it may remain trapped (Schorghofer & Taylor, 2007). The Oxford 3D Thermophysical Model (Warren et al., 2019) was applied to the PSRs studied in this work to estimate the temperature of the regolith with depth. The model analyzes a grid at each location with 8 ppd in the longitude direction and 16 ppd in the latitude direction. If the maximum temperature below the surface remains below 112 K, it is assumed that stable water ice could be found. The results show that water ice can be found at depths below 0.05 m in one PSR, and below 0.10 m in three further PSRs in this study (see Figure 5 and Table S6). The remaining PSR in this study was too small for the Oxford 3D Thermophysical model to be applied without the illuminated regions being included in the analysis. The estimated ice stability depths combined with the bearing capacity data suggest that low-latitude PSRs that have the potential to contain buried water are likely traversable. Further work will be needed to assess bearing capacities and, thus, trafficability in the higher latitude, colder PSRs. In situ analyses are also needed to confirm whether water ice is present at such locations and how it is distributed.

Technique Limitations
It is important to understand the limitations of the measurements and the application of the Hansen (1970) bearing capacity equation, which were discussed in detail in Bickel et al. (2019). A number of assumptions were made when determining the effective incidence angle, θ, of secondary illumination of PSRs, which affects estimated track depths and bearing capacity (see Figure S3). The uncertainty in θ is not well defined, but to illustrate the consequences of that uncertainty, an error analysis was performed assuming 5°and 10°e rrors in θ. For depths <6 m, a 5°and 10°uncertainty in θ results in a maximum uncertainty in our calculated bearing capacity values of q þ107 f −95 kN/m 2 and q þ59 f −51 kN/m 2 , respectively. The largest impact of an uncertainty in θ occurs when measurements were taken at low slope angles, where the resulting θ is small. This is a consequence of the trigonometric nature of the θ estimations (see Figure 4c).
Two additional sources of potential uncertainty lie with the track depth estimations. (1) When estimating track depth, a first-order approximation was made, which assumes a flat-bottomed track. However, boulder tracks are likely to be irregular in shape, particularly as lunar boulders are particularly angular Kumar et al., 2016). As a consequence, the track depth estimations may be shallower than the actual tracks. (2) At very high incidence angles, the rim of boulder tracks might cast a very long shadow that covers a majority of the track. These long shadows do not allow for an estimate of depth at the center of the track as desired, but toward the opposite rim, ultimately resulting in inaccurate track depth estimates. The consequence is a bearing capacity that is associated with a track on average 23% deeper than originally thought.

Implications for PSR Regolith Processes
Boulder track depths measured in this work range from 0.28 to 4.68 m, and, therefore, the estimated bearing capacity results are only directly applicable in that range. When comparing the measurements reported here in PSRs with those in mare and highland regions , the bearing capacities are statistically the same in all regions in the 1-to 2-m depth, which indicates that regolith processes at this depth are unaffected by temperature differences at the surface. The bearing capacities are statistically higher for PSRs at the shallowest depths (≤1 m) as compared to highland and mare regions. Interestingly, the bearing capacities in the PSRs at ≤1 m are similar to those in pyroclastic deposits at equivalent 10.1029/2019JE006157 Journal of Geophysical Research: Planets depths. One might posit that extra strength is due to surficial or cohesive intergrain ice present in PSRs. Although remote sensing evidence suggests that there is no surficial water ice, there could be ice at depths of >10 cm as suggested by the Oxford 3D Thermophysical model. Water could have migrated beneath the surface as a consequence of impact bombardment and regolith overturn (Hurley et al., 2012). Providing water has migrated deeply enough into the regolith, perhaps the presence of water ice connects the grains forming a rigid matrix structure. If the PSRs studied in this work do not contain water ice, then an explanation is required for the statistically higher bearing capacity of PSR regolith in the upper meter. The specific location of measurements may be a factor in the increased PSR regolith strength calculated using this method. It is presumed that PSR regolith is porous; however, this study considers boulders and tracks, which are generally located on a slope. Slow downslope movement of regolith over time (see Bickel et al., 2019;Carrier et al., 1991) might have resulted in a modified porosity and density of the regolith on the slopes. These modified properties might vary from those at or closer to the center of PSRs. In situ measurements would be required to understand how representative the bearing capacity results are as compared to other locations within and around PSRs at different latitudes.
If subsurface ice existed at relatively shallow depths when the tracks were produced, the ice may have ablated from the walls of the tracks. Based on studies of boulder tracks elsewhere (Hurwitz & Kring, 2016;Kumar et al., 2016), tracks are obliterated by regolith processes on time scales of 20 to 35 Ma, although icy regolith may result in different track aging rates. Experimental data suggest that ice sublimates from hydrated lunar regolith simulants at a rate of ≤2% in 5 days at temperatures of~100 K (Piquette et al., 2017). Therefore, any ice exposed by the rolling boulders may have ablated before the tracks were observed by LROC.
The dimensions of the boulders measured in this work range from 4.2 to 16.3 m. This is a significantly smaller range than that found in highland, mare, and pyroclastic regions , which contain boulders of 1 to 31 m in diameter (see Figure S2). The lack of small boulders in PSRs may be an artifact of image noise, which makes it difficult to identify small features. The reason for a smaller proportion of larger boulders, relative to that seen in mare and highland regions, is less clear.

Implications for PSR Traverses
The observations indicate that rovers with wheel radii of ≳~28 cm should be able to traverse PSRs of the type studied here. It remains possible that the upper few centimeters of regolith are physically different than deeper regolith analyzed in this study. For example, the upper few centimeters of lunar regolith are generally very porous at~52% (Carrier et al., 1991); however, if water ice is present, then the porosity could be decreased if the ice fills the pore spaces. Additional work may, therefore, be needed to determine if vehicles with wheel diameters ≲28 cm can traverse these regions. The comparative assessment of average PSR bearing capacities with those of other terranes used highland and mare measurements obtained mostly at high slope angles (>15°). Because slope affects bearing capacity, that may have produced a bias toward lower bearing capacity values in highland and mare terranes. That prompts a second comparison of PSR properties with those of highland and mare using only those observations made on slopes >15°(see Table S5). On those surfaces, PSR regolith appears to be statistically stronger in the ≤1-m range, with the caveat that the result is based on only two data points. Both comparative approaches have weaknesses but point to similar conclusions: that PSRs appear to be stronger than highland and mare regions at equivalent depths and slope angles. When comparing the results of PSRs with pyroclastic deposits, pyroclastic deposits appear to be statistically similar in strength to PSRs when all slope angles are considered, but statistically weaker when considering higher slope angles only. As the pyroclastic data have a less variable distribution across all slope angles, it is more appropriate to consider the statistical analysis for all slope angles, which indicates that pyroclastic deposits are similar in strength to PSRs at ≤1 m. Despite those trends, it is also true that there is a range of bearing capacities in any terrane.

Implications of Regolith Relative Density on Bearing Capacity
The value for bulk density of PSR regolith used in this work (1.66 g/cm 3 for regolith >30 cm deep) was selected because the lunar polar regolith is thought to be mostly feldspathic "highland-like" material. Therefore, the bulk density value is taken from , which analyzed highland material. However, it has also been theorized that PSRs are highly porous , which 10.1029/2019JE006157 Journal of Geophysical Research: Planets would equate to a lower relative and bulk density. Bulk density, ρ (g/cm 3 ), and porosity, n, are related as follows (Carrier et al., 1991): where G is the specific gravity of regolith taken to be 3.1 and ρ w is the density of water (1 g/cm 3 ). When using a bulk density of 1.66 g/cm 3 , the estimated porosity would therefore be~46%. As evidence suggests that PSRs could have porosities as high as 70% , the equivalent bulk density is calculated to be 0.93 g/cm 3 . As illustrated by Carrier et al. (1991) using a basaltic regolith simulant, an increase in porosity and a decrease in relative and bulk density will result in reduced cohesion and friction angle values of the respective soil. However, as the input parameters for cohesion and friction angle adopted from Bickel et al. (2019) are already conservative assumptions, particularly the angle of internal friction (29°), a further reduction would likely result in unrealistic bearing capacity estimates. A comparison of the bearing capacities as calculated using this lower-limit value for bulk density is shown in Figure S4 and Tables S7 and S8.
The results show that the bearing capacity is reduced at higher porosities (e.g., from 127 ± 29 to 73 ± 14 kN/m 2 in PSRs at depths of ≤1 m and on a 0°slope for densities of 1.66 and 0.93 g/cm 3 , respectively); however, the upper meter of regolith in PSRs is still estimated to be as strong as highland regolith and equally as strong as mare and pyroclastic regolith. At depths >1 m, the bearing capacity of regolith in PSRs is calculated to be significantly lower (~34-45% lower) than highland, mare, and pyroclastic regolith for the lowest bulk density scenario; for example, PSRs have a bearing capacity of 110 ± 16 kN/m 2 when applying a density of 0.93 g/cm 3 as compared to bearing capacity values calculated for highland, mare, and pyroclastic regions of 166 ± 24, 167 ± 19, and 200 ± 30 kN/m 2 , respectively, at depths of ≤1 m and on a 0°slope. As a consequence, the PSR regolith may be somewhat weaker at depths of >1 m, as compared to other regions of the Moon.

Requirements for Increasing the PSR Boulder Data Set
Although the data support a self-consistent set of conclusions, we recognize that it needs to be confirmed. Now that the methodology has been demonstrated, the data set could be expanded. In this work, 13 boulders across five locations were analyzed. It seems reasonable that more boulders and boulder tracks can be located and measured within the north polar region in a similar latitude range of~70-75°. That will allow an assessment of heterogeneities that may exist in PSRs of different ages, composed of different lithologies. Some progress should be possible with the existing set of LROC NAC images, but we envision that a number of boulder tracks will be captured in some of the darkest and coldest PSRs once ShadowCam is deployed (Robinson & ShadowCam Team, 2018). However, as the spatial resolution of ShadowCam will be lower than the LROC imager, it will likely only be able to resolve large boulder tracks, which cannot be used to estimate bearing capacity at shallow depths.

Conclusions
A method was developed to study boulder tracks in NAC images of lunar south polar PSRs from which the bearing capacities of PSR regolith could be assessed qualitatively and quantitatively. Boulders that traveled from sunlit regions into PSRs are seemingly unaffected. They do not disappear into an overly porous or weak regolith. Rather, the dimensions of the boulder tracks remain unchanged suggesting that there are no large bearing capacity contrasts between the regolith in sunlit and permanently shaded regions. Those track dimensions can also be used to estimate bearing capacities. When compared with bearing capacities determined in other terranes , regolith in PSRs appear to be as strong as the regolith in highland and mare regions, and similar to that in pyroclastic deposits within the upper~0.28 to 1 m of regolith across all slope angles. Those qualitative and quantitative results seem to contradict reports of very high porosity in PSRs (e.g., Gladstone et al., 2012;Hayne et al., 2017;Metzger et al., 2018;Schultz et al., 2010) or require those proposed high-porosity conditions be limited to the uppermost 28 cm or the presence of water ice provides additional strength in a porous regolith. In either case, trafficability of PSRs may be possible with rovers of wheel radii ≳28 cm with sufficient traction. In situ analyses are still required to verify these results and to establish the distribution of regolith strength within the uppermost centimeters and meters, and the lateral heterogeneity of the regolith, of PSRs.