Small-scale topographic irregularities on Phobos: image and numerical analyses for MMX mission

The mothership of the Martian Moons eXploration (MMX) will perform the first landing and sampling on the surface of Phobos. For the safe landing, the 2.1-m-wide mothership of the MMX should find a smooth surface with at most 40 cm topographic irregularity, however, whose abundance or even existence is not guaranteed based on current knowledge. We studied the highest resolution (a few meters per pixel) images of Phobos for possible topographic irregularities in terms of boulder (positive relief feature) and crater distributions. We find that the spatial number densities of positive relief features and craters can vary significantly, indicating that the surface irregularities vary significantly over the entire surface. We extrapolate the size-frequency distributions of positive relief features to evaluate the surface roughness below the image resolution limit. We find that the probabilities that topographic irregularities are  < 40 cm for the areas of 4 × 4 m and 20 × 20 m are  > 33% and  < 1% for boulder-rich areas and  > 88% and  > 13% for boulder-poor areas, respectively, even for the worst-case estimates. The estimated probabilities largely increase when we reduce the assumed number of positive relief features, which are more realistic cases. These indicate high probabilities of finding a smooth enough place to land on Phobos’ surface safely.


Introduction
Phobos and Deimos, the two moons of Mars, are intriguing not only for their high scientific interests, but also as targets of potential human missions. Martian Moons eXploration (MMX) is the Japan Aerospace Exploration Agency (JAXA)'s mission to explore Phobos and Deimos, scheduled to be launched in 2024. The MMX spacecraft will perform in situ observations of both Phobos and Deimos, land and collect samples on Phobos, and bring them back to Earth. In addition, the MMX spacecraft will deploy a small rover developed by Centre National d'Etudes Spatiales (CNES) and the German Aerospace Center (DLR) for in situ investigations of the surface properties.
Designs of the mothership, the rover, and the sampling device depend largely on the target body's surface conditions. Thus, the Landing Operation Working Team (LOWT) and the Surface Science and Geology Sub Science Team (SSG-SST) of MMX are organized by scientists and engineers to evaluate Phobos' surface conditions. The team's views for both the surface environments and regolith properties are summarized in the accompanying paper .
Previous Phobos' explorations provide useful information for the MMX mission, including photomosaics, maps, and numerical shape models (Basilevsky et al. 2014;Karachevtseva et al. 2014;Oberst et al. 2014;Salamunićcar et al. 2014, Wählisch et al. 2014Willner et al. 2014). However, to evaluate the landing hazard, we need to evaluate the topographic irregularities at a scale of the order of tens of centimeters, which is difficult to do based only on the current numerical shape model or available images at a lower resolution than needed. Thus, we attempt pixel-scale image mapping and numerical evaluations based on the mapping results to evaluate the worst-case topographic irregularities.
In this work, we provide a brief overview of technical characteristics and restrictions of the MMX mission to Phobos' surface ("Topographic irregularities and the engineering safety" Section) and summarize typical geological features of Phobos (i.e., topographic obstacles) identified by previous studies ("Previous observations" Section). To estimate surface roughness at a high spatial resolution, we perform mappings of craters and boulders at a sub-pixel scale ("Analyses of high-resolution images" Section) and artificially develop numerical terrain models by extrapolating the estimated boulder size-frequency distributions ("Artificial terrain model" Section). Then we measure the topographic irregularities of these terrain models to assess landing safety ("Results and discussions" Section). We conclude that there are plenty of smooth enough areas favorable for landing on Phobos ("Conclusion" Section).

Topographic irregularities and the engineering safety
In this paper, the term "topography" refers to the surface forms and features themselves. We use the topographic height for the topographic difference (or approximately the elevation difference) for a topographic feature's normal distance from the reference plane defined in an area. Note that elevation is usually defined as the height above or below a reference point, mostly a reference geoid, which is difficult to define for Phobos under the current understanding of its internal densities.
Phobos is a small satellite of 26 × 22.8 × 18.2 km in size. Like a similarly sized asteroid Eros, Phobos' major landforms are craters and grooves with many boulders distributed over the body. Most distinctive is the Stickney crater, which is about 8 km in diameter. It is similar in size to the Shoemaker crater on Eros.
Based on spacecraft's observational data, including Mariner-9, Viking-Orbiter, Mars Global Surveyor, Mars Reconnaissance Orbiter, and Mars Express, numerical shape models have been built by several researchers leading to the knowledge of the overall shape of Phobos (e.g., Willner et al. 2014). We note that, even with a numerical shape model, determining the local gravity direction is not simple due to the Phobos' proximity to Mars. Nevertheless, we can calculate the gravity field with a constant density value for Phobos' entire body to evaluate the local gravitational acceleration in every facet (i.e., the face that constitutes the polygon shape model). We find that the slopes with respect to the local gravities are mostly lower than 40° (Wang and Wu 2020;Willner et al. 2014), with plenty of < 10° slope areas , which are preferred for a lander and a rover. These values are useful for designing operational schemes, landing procedures, and selecting preliminary landing sites.
The spacecraft's size is about 2.1 m in width (excluding the solar paddles), and the vertical distance between the tip of the probe's legs and the mechanism carrying the various instruments is about 40 cm. Thus, acceptable topographic irregularity for its safely landing is only 40 cm at least under the current design. Considering the extension of the landing pads from the spacecraft, we needed to evaluate the surface roughness for areas of about 4 × 4 m as possible landing areas. The original operation plan for the landing of the MMX mothership is as follows: (1) descending the spacecraft's altitude from a Quasi-Satellite Orbit (Kuramoto 2021) to about 1 km in altitude to eliminate the horizontal velocity with respect to the surface; (2) descending the altitude by using onboard terrain matching system to the altitude where a target-marker is released; (3) further lowering the spacecraft's altitude to 20 m or so; and (4) performing an almost free-fall touchdown without thruster firing to prevent plume contamination. The accuracy of the horizontal drifting suppression is one of the designing parameters, which significantly controls the landing accuracy. Thus, topographic irregularities in varying areal scales, such as 10 × 10 m-40 × 40 m, are needed to be evaluated.
The averaged facet area of the available shape model with the highest resolution is about 10 4 m 2 , which exceeds the scale discussed above. Also, even with the highest resolution image of Phobos, which is about a couple of meters per pixel, evaluations of topographic irregularities at the required high resolutions are difficult because at least a few to ten pixels are needed to recognize a topographic feature. Thus, we need further efforts to assess the surface roughness at a higher resolution for landing hazard evaluation.

Previous observations
The surface of Phobos was globally imaged by high-resolution stereo camera (HRSC) of Mars Express, resulting in a global mosaic with a resolution of about 12 m/ pixel (Wählisch et al. 2010(Wählisch et al. , 2014. The global image (e.g., Fig. 1a) provides essential information on the global distribution of geological features. Detailed surface features are also studied in some limited areas. The super-resolution technique allows a resolution of the selected area of HRSC to be about 2 m per pixel. The highest resolutions achieved on images of Mars Orbiter camera (MOC) aboard Mars Global Surveyor are about 1-7 m per pixel for a few areas.
Based on these images, geological studies find that dominant topographic landforms on Phobos' surface differ in scale and location and are typically classified into three major morphologies: craters, grooves, and boulders (Basilevsky et al. 2014). Craters are commonly distributed over the body, and their diameters range from at least several meters to ~ 10 km (Salamunićcar et al. 2014;Schmedemann et al. 2014;Hartmann and Neukum 2001).
A total of 1072 craters larger than 250 m in diameter were identified on Phobos , whose statistics are close to the Moon's highlands (Thomas and Veverka 1980). Some craters appear to be fresh, while some are highly degraded. Most of their depth-to-diameter ratios are between ~ 0.02 and ~ 0.2 (Hemmi and Miyamoto 2020;Basilevsky et al. 2014;Karachevtseva et al. 2014).
The smallest craters (about 3 m in diameter) were mapped on the highest resolution image (SP2-55103) obtained by the MOC. Considering that much smaller craters down to 0.1 μm (Hörz et al. 1975) exist on the surface of the Moon, which are not visible in the images of Phobos, cratering events down to micrometers should be present on the surface of Phobos (Basilevsky et al. 2014). Crater equilibrium was suggested for certain diameter ranges (e.g., Hartmann and Neukum 2001;Thomas and Veverka 1980). However, the actual spatial distribution and size-frequency distribution (SFD) of craters smaller than a few meters in diameter have not been quantitatively characterized due to the image resolution limit.
Grooves and pit chains exist as linear depressions, several of which are parallel to each other. Groove sizes vary from 23 to 475 m in width and from 2 to ~ 30 km in length (Murray and Heggie 2014). Nearly five hundred grooves are distributed globally, except for the trailing hemisphere of Phobos (e.g., Kikuchi and Miyamoto 2014;Thomas et al. 1979;Murray et al. 1994). Their nature is of high scientific interest, and their origins may hold important keys in the evolutional and structural history of Phobos. However, the landing operation team will most likely avoid landing the spacecraft into the bottom of a groove due to engineering safety reasons. Thus, we neglect the influence of grooves on the surface irregularities in this paper.
Numerous boulders are identified in high-resolution images of Phobos. They range from tens of meters down to a few meters in diameter. The nature of boulders smaller than a few meters in size is uncertain due to the paucity of higher-resolution images taken at the moderate incidence and phase angles (Thomas et al. 2000;Karachevtseva et al. 2014).
Previous geological studies of meter-scale features are limited (Basilevsky et al. 2014(Basilevsky et al. , 2015Karachevtseva et al. 2014) due to the limited availability of high-resolution images. A high-resolution image of an area at least wider than about 1 km 2 is needed for the statistical study of the distributions of meter-scale craters and boulders. Considering the spatial resolution, the spatial extent of illuminated areas, and the degree of noise and blurs, among the available datasets of Phobos, the Mars Global Surveyor MOC image SP2-55103 is the most suitable for the analyses of numerous boulders and craters as studied by Karachevtseva et al. (2014). Thus, first, we perform a similar but more detailed study of the same image to obtain statistical information.

Analyses of high-resolution images
We mapped both craters and positive relief features down to sub-pixel scales on MOC image SP2-55103, which is the highest-resolution image of the sub-Mars surface of Phobos of an area wider than several km 2 . The ground pixel scale of this image is about ~ 2.43 m on average, but we focus on regions with the highest resolution of about 1.1 m (due to the irregular shape of Phobos in this image frame, the image resolution largely varies even within the image). Boulders smaller than several meters are difficult to identify even in the highest resolution images. Thus, we picked up all positive relief as boulder candidates, which are composed of bright pixels (in the sunny side) immediately next to dark pixels (in the shadow side) along the line of the solar azimuth angle.
We selected two regions of study (Regions A and B) for precise mapping (Fig. 1b). Such regions are selected because Region A has relatively higher densities of positive relief features and craters, while Region B has relatively lower densities of those same features than surrounding regions. Note that Region B is composed of 4 subregions (Regions B1-B4) to avoid the shadowed area inside a large crater within somehow similar illumination conditions as Region A. Regions A and B (the sum of Regions B1 to B4) have comparable areas, ~ 1.14 km 2 and ~ 1.20 km 2 , respectively (the individual areas of Regions B1-B4 are within ~ 0.26-0.31 km 2 ). We measured diameters and center coordinates of craters and positive relief features in these regions by using the Cra-terTools v2.1 extension (Kneissl et al. 2011) for the ESRI ArcMap 10.1 or higher. To minimize the distortions of measurement areas due to map projection of irregularshaped Phobos' surface (deviated from an 11.1 km-radius spherical reference), we applied the CraterTools' Topography Correction (Kneissl et al. 2014) with the use of the global digital terrain model of Willner et al. (2014) to our mapping results. The mapped craters and positive relief features are categorized into either "confirmed" or "candidate", depending on their confidence levels. Confirmed features are those with evident morphological characteristics identified from multiple solar-incidence images [HiRISE image PSP_007769_9015 and the Stooke map ( Fig. 1a)] by several researchers. However, candidate features are either nearly sub-pixel-scale or highly degraded features, recognized only from patterns of positive/negative relief ( Fig. 2a-d; e.g., the arrangement of darker, brighter, then much darker pixels along a subsolar azimuth direction implies the presence of one or subpixel-scale positive relief feature at the brighter pixel).
We study the local incidence angles of the MOC image, derived from the bundle-adjusted MOC camera position/pointing information and the shape model of Phobos with 100 m/pixel resolution ). The incidence angles range from ~ 70 to ~ 80 degrees within the studied regions (Fig. 1c, d). This greatly helps to identify relatively low topographic features. Thus, although candidate features may include imaging artifacts/noises, most would represent actual positive relief features. In this sense, we assume that the mapping result is worstcase. Also, the difference in the mean solar-incidence values between Regions A and B is only about 4 degrees (Fig. 1e), which is expected to provide a very similar or slightly better appearance of topographic features in Region A than Region B (e.g., a 2-m high feature will give ~ 6.2 and ~ 8.0-m-long shadows at solar incidence angles of 72 and 76 degrees, respectively, which result in at most one-pixel-wide difference between two regions in the MOC image).
We find that the spatial number densities of craters and boulders can vary significantly even in the same image frame (Fig. 3); Region A has 202 km −2 confirmed positive relief features (and 2277 km −2 including candidates) and 278 km −2 confirmed craters (1226 km −2 including candidates), while Region B has 46 km −2 positive relief features (224 km −2 including candidates) and 234 km −2 craters (460 km −2 including candidates). The mapping results of both boulders and craters are statistically summarized as cumulative plots in Figs. 4, 5. We define their size-frequency distributions (SFDs) following the conventional representation, the total number per unit area (N) of target features (boulders or craters) larger than diameter D, i.e., N = N D0 (D⁄D 0 ) −α , where D 0 is a reference diameter and N D0 is the N value of boulders/craters with a diameter larger than D 0 (10 m for boulder SFDs and 1 km for crater SFDs). For craters, a power-law exponent α = 2 is generally considered to represent an equilibrium population (Melosh 1989).
The SFDs of confirmed positive relief features are shown as red dots in Fig. 4a, b in a log-log plot (the slope of the distribution lines in such a plot corresponds to the power-law exponent in the actual distribution). Note that, in the range of ~ 5-8 m in diameter,  the cumulative slope (fitted to about − 6) is similar to that of Karachevtseva et al. (2014). The SFDs of both confirmed and candidate positive relief features (yellow dots) appear to follow the slope of ~ − 3.2 in the range of 1-2 m in diameter. The SFDs of confirmed craters and the sum of both confirmed and candidate craters are shown in Fig. 5. In these log-log plots, the SFDs roughly follow the lines with a slope of − 2 at D = ~ 10-30 m, which are consistent with the plots of the distributions of larger diameter craters measured by previous studies (Salamunićcar et al. 2014;Hartmann and Neukum 2001). The equilibrium crater density (a few to 10 percent of geometric saturation) (Gault 1970) seems to be achieved.

Artificial terrain model
We have identified the locations and diameters of confirmed craters, crater candidates, confirmed positive relief features, and positive relief feature candidates through the analysis discussed above. Assuming that the cm to meter-scale topography on Phobos mostly results from craters and boulders distributions, we might evaluate the topographic roughness by creating an artificial terrain model (DTM) based on such information.
To artificially develop a DTM, we define two-dimensional grid arrays of 10,000 × 10,000 elements for simulating a region of 1 × 1 km, which means the resolution of this model is 10 cm/element. Each element has its topographic height value, whose initial value is 0.0 m (i.e., a reference level). We assume that each crater is composed of a circular uplifted rim around a bowlshaped depression and their morphology corresponds to the one given by a standard crater gravity scaling law. We put such craters in the terrain model by subtracting the height values corresponding to the area of the bowl-shaped depression and adding for the circular uplifted rim. We place the craters in descending order; larger craters are placed earlier.

Diameter [m]
10 0 10 1 10 2 10 3 10 -2 10 -1  By using Geospatial Data Abstraction Library (GDAL) version 3.1.2 and the ArcGIS's Spatial Analyst tool, the resulting artificial DTMs were converted to shaded relief images with shadows, which are artificially illuminated at an azimuth of 270° (from the left of the figure to the right) and an altitude of 15° (i.e., the roughly similar illumination condition of the target region (Fig. 6a) in the MOC image SP2-55103). Figure 6b shows the shaded relief image of the resulting artificial DTM with confirmed craters, whose d/Ds are randomly given between 0.01-0.2 according to the simple uniform distribution. Note that, even though we used both crater locations and diameters from mapping results, the resulting DTM is not similar to the original MOC image. This result indicates the importance of evaluating the d/D (depth to diameter) of craters.

Crater depth-to-diameter ratio
The topographic characteristics of craters below 100 m in diameter are difficult to evaluate from the numerical shape (e.g., Willner et al. 2014). Basilevsky et al. (2014) partially overcame this issue by carefully evaluating images to obtain craters' 2D profiles. They classified craters into three morphologic classes: Class 1 for those with > 0.1 in d/D (depth to diameter) with steepest inner slopes (> 20°), Class 2 for those with 0.05-0.1 in d/D with shallower inner slopes (10-20°), and Class 3 for those with < 0.05 in d/D with shallow inner slopes (< 10°). Such variations may be the results of degradations through many possible processes.
Following this idea, we assume that each crater initially has a traditional simple bowl shape. The simple crater's rim-to-floor depth is about 1/5 of its diameter, and its sharp-crest rim stands about 4% of the crater diameter above the surrounding crater (Richardson 2009). The change in elevation for the crater interior and exterior is where r is the radial range from the crater's center, R is the rim crest radius of the crater, h r is the rim height, d is the depth of the crater, and η is the radius of the continuous ejecta blanket (in our model, η = 4 as the continuous ejecta blanket extends to 4 crater radii). A degradation process is simulated with the following two-dimensional diffusion equation (e.g., Richardson et al. 2020): where κ is diffusivity. The degradation process is simulated by fixing the coefficient and setting the appropriate diffusion time t. As a result, we obtain 20 profiles of the crater's morphology with its d/D of 0.01-0.2. The spline interpolation method is used to keep the continuous shape of the resulting crater morphology (Fig. 7).
We modified the d/D of the confirmed craters to make the appearance of DTM similar to the MOC image. Note that this procedure is similar to the shapefrom-shading method in some sense, but is also different because we empirically assume the shapes of craters and neglect any possible brightness difference due to surface textures. Thus, we can focus on each shadow length of each crater to adjust d/D. For a deep crater with its shadow cast from its wall to its interior, we selected the profile (d/D) of a crater that shows the same shadow length as the one observed in the MOC image. For a shallow crater without any interior shadow, we used the profile of the deepest shadowless crater at a given solar incidence angle at its location because the worst-case crater should be considered for the assessment of future landing operation. Figure 6c shows the result of DTM with modified d/D for all confirmed craters. As for crater candidates, d/D should generally be very small because otherwise we can identify them as a confirmed crater. Thus, we chose 0.01-0.02 for d/D of the candidate craters and placed the craters in Fig. 6c. The shaded relief image of the resultant DTM is shown in Fig. 6e.
Note that we excluded craters smaller than 82.7 cm in diameter in our artificial models. This comes from the aim of this study, i.e., the assessment of landing safety, most of which is based on Rodgers et al. (2016). For the landing performance of the lander footpads, one small boulder is more hazardous than one small crater. A footpad placed on a single acuate (convex-upward) boulder with a certain height is very unstable; however, a footpad over/inside a small (1-2 × footpad-wide) crater is in a (quasi-)stable state because of the convex-downward geometry of a crater profile. Besides, little is known about the spatial/size-frequency distributions of submeter craters on small bodies, compared to those of sub-meter boulders on small bodies (e.g., asteroid Eros, Itokawa, Ryugu, etc.). For the reasons above mentioned, we focused on the distributions of positive relief features only for the sub-meter-scale elements of our artificial DTMs.

Shapes of positive relief features
We put positive relief features on the DTM as bowlshaped uplifts. Similar to craters, shapes of positive relief features are difficult to define. For simplicity, we use an ellipsoid to represent a positive relief feature. Thomas et al. (2000) measured the boulders' width/height ratio by measuring the heights from lengths of shadows, which showed that the height/width ratios follow a normal distribution with its mean 0.25 and standard devia-

Extrapolation of positive relief feature's SFD
The resultant DTM (Fig. 6f ) shows a consistent appearance to Region A's original MOC image. Because the illumination condition of Fig. 6e is carefully arranged to be the same as that of the original image, the overall similarities may indicate that we can evaluate surface irregularities if we can appropriately estimate the distributions of smaller blocks below the resolution image. However, extrapolation of an SFD is very challenging because the observed SFDs on small bodies appear to saturate at different sizes, which may reflect some nature of granular materials. However, these could be biased by resolution limits. Thus, to make the discussion simple, we simply extrapolated our SFDs from the analysis discussed above, knowing that it overestimates the numbers of smaller particles. We consider this is still justified because (1) our primary purpose is for engineering safety evaluations, and (2) the number of particles is always overestimated, which always gives the worst case.
We prepared two types of model terrains, Models 1 and 2, which correspond with two target regions, Regions A and B, respectively. Small positive relief features, down to 35 cm, were placed spatially randomly. The extrapolations of their SFDs of Models 1 and 2 were made using

Results and discussion
We artificially develop DTMs of Region A in Fig. 1 by putting confirmed craters, crater candidates, confirmed positive relief features, positive relief features candidates, and small positive relief features. The shaded relief images of the DTMs developed in this study are shown in Fig. 8c, d.
We study the frequency of the maximum topographic differences within a given small region as shown in Fig. 9a, d for Models 1 and 2, respectively. Naturally, the larger the size of the given region, the higher the peak of the maximum topographic difference. The cumulative probabilities of the maximum topographic difference are also calculated (Fig. 9c, d).
We find that the probabilities that topographic irregularities are < 40 cm for the areas of 4 × 4 m and 20 × 20 m are > 33% and < 1% for Model 1, and > 77% and > 2% for Model 2. This indicates that, even in the worst case, we can still find a smooth area within 40 cm in topographic differences with a probability of at least 30%. Thus, even if we do not control the spacecraft at the time of landing, the probability of successful landing in a smooth region is more than 30%.
As we discussed above, we adjusted the d/D of craters to simulate the original MOC image appropriately. The resultant d/D distribution is shown in Fig. 10. We find that this is a useful set to develop DTMs artificially. We use the statistical distribution of d/D in this figure, we randomly plot confirmed craters and crater candidates on an initially flat DTM to simulate Region B, which is a boulder-poor region. The number of positive relief features is also extrapolated down to 35 cm in diameter as for Region A as shown in Fig. 11a, b. The shaded relief images of the resulted DTMs are shown in Fig. 11c, d. Their statistics are shown in Fig. 12. We find that, in this case, the probabilities that topographic irregularities are < 40 cm for the areas of 4 × 4 m and 20 × 20 m are > 79% and > 2%, respectively, for Model 1 and > 88% and > 13%, respectively, for Model 2.
We note that these two areas could be considered as boulder-rich regions on the Phobos' surface. We surveyed other areas far from Regions A and B using relatively high-resolution images taken by the HRSC on the Mars Express. We find that, in some regions such as those in Fig. 13, however, only a few boulders were identified. This indicates that the areas other than the eastern side of Stickney crater may be more favorable in terms of topographic irregularities for landing.

Conclusion
Evaluation of topographic irregularities below image resolution is very challenging, although it is required for the appropriate design of the landing of the MMX spacecraft. Thus, we study the topographic irregularities of Phobos' surface based on the area where the highest resolution image is available. By combining image analysis and numerical development of DTMs at high resolution, we find that Phobos' surface likely has plenty of smooth areas that are favorable to the landing of the MMX spacecraft. However, the maximum topographic difference increases with the area width, which means that the landing safety requires a landing accuracy that is optimized for the estimated width of smooth areas.