Blocks Size Frequency Distribution in the Enceladus Tiger Stripes Area: Implications on Their Formative Processes

We study the size frequency distribution of the blocks located in the deeply fractured, geologically active Enceladus South Polar Terrain with the aim to suggest their formative mechanisms. Through the Cassini ISS images, we identify ~17,000 blocks with sizes ranging from ~25 m to 366 m, and located at different distances from the Damascus, Baghdad and Cairo Sulci. On all counts and for both Damascus and Baghdad cases, the power-law fitting curve has an index that is similar to the one obtained on the deeply fractured, actively sublimating Hathor cliff on comet 67P/ChuryumovGerasimenko, where several non-dislodged blocks are observed. This suggests that as for 67P, sublimation and surface stresses favor similar fractures development in the Enceladus icy matrix, hence resulting in comparable block disaggregation. A steeper power-law index for Cairo counts may suggest a higher degree of fragmentation, which could be the result of localized, stronger tectonic disruption of lithospheric ice. Eventually, we show that the smallest blocks identified are located from tens of m to 20–25 km from the Sulci fissures, while the largest blocks are found closer to the tiger stripes. This result supports the ejection hypothesis mechanism as the possible source of blocks.


Introduction
Enceladus is a~500 km-size icy moon of Saturn [1], orbiting at a distance of 3.94 Saturn radii from the planet. Its surface is dominated by pure water ice [2], heavily cratered in the northern latitudes as well as in its trailing and leading mid-latitude hemispheres. Through Voyager 1 and 2 km-scale images [3,4], it was assessed that Enceladus is characterized by different geological units, documenting a complex origin and evolution.
In 2005, the Cassini Imaging Science Subsystem-Narrow Angle Camera (ISS-NAC, [5]) acquired high-resolution images of the moon's surface [6], revealing an almost craterless, fracture-dominated and geologically active province surrounding the South Pole of Enceladus, called the South Polar Terrain (SPT). Inside this unit (called csp, i.e., central south polar, surrounded by the southern curvilinear unit, called cl3, [7]), both cryovolcanic and tectonic activities are evident [7][8][9], with icy geysers emanating water vapor plumes from four, 100 km-long equally-spaced tension fractures [6], called "tiger stripes" (hereafter named TS), Figure 1.  [7]. The red rectangle shows the location of our study area, presented in Figure 2.  [7]. The red rectangle shows the location of our study area, presented in Figure 2.  [7]). The yellow outline shows the full area covered by the dataset. On the left side of the saw-tooth yellow boundary, a bad signal-to-noise part of the image has been eliminated. The red rectangle shows the location of Figure 3, while the orange circles show the jets' sources identified by Porco et al. [10]. Porco et al. [10] evidenced that the origin of such jets is strictly related to normal tensile stresses that open vertical pathways inside the TS (the fractures open and close from apocenter to pericenter, [11]), reaching the ocean that is located between 5 to 40 km  [7]). The yellow outline shows the full area covered by the dataset. On the left side of the saw-tooth yellow boundary, a bad signal-to-noise part of the image has been eliminated. The red rectangle shows the location of Figure 3, while the orange circles show the jets' sources identified by Porco et al. [10].  [7]). The yellow outline shows the full area covered by the dataset. On the left side of the saw-tooth yellow boundary, a bad signal-to-noise part of the image has been eliminated. The red rectangle shows the location of Figure 3, while the orange circles show the jets' sources identified by Porco et al. [10]. Porco et al. [10] evidenced that the origin of such jets is strictly related to normal tensile stresses that open vertical pathways inside the TS (the fractures open and close from apocenter to pericenter, [11]), reaching the ocean that is located between 5 to 40 km Porco et al. [10] evidenced that the origin of such jets is strictly related to normal tensile stresses that open vertical pathways inside the TS (the fractures open and close from apocenter to pericenter, [11]), reaching the ocean that is located between 5 to 40 km underneath [12][13][14][15]. Through these conduits water vapor and liquid can reach the surface, hence depositing the latent heat they have [16,17].
In order to have higher-resolution views of the geyser vents and to study the TS rugged floors and their closest surroundings, three Enceladus close flybys have been performed by Cassini in the August-October 2008 timeframe, and in November 2009. The returned m-scale images have been unprecedented and led to the identification that the spacing of geysers is uniform along the three most active tiger stripes (Damascus, Bagdad and Cairo Sulci, [18]). In addition, the orientation of the observed jets is strongly correlated with the directions of the TS, crosscutting fractures or the orientation of local tectonic features [18].
The meter-to-decameter scale ISS-NAC images revealed also the presence of positiverelief icy blocks located in the SPT (hereafter called blocks), situated in close proximity to the TS. Martens et al. [1] focused on the spatial density of such features, showing that the variations in number density can be substantial over short spatial scale. Since impact cratering cannot be considered as the source of the SPT blocks, nor seismic disturbance or mass wasting, Martens et al. [1] suggested that sublimation processes, ballistic emplacement through geysers eruption and/or tectonic disruption of lithospheric ice can be the source of the blocks.
Besides the blocks spatial density, it is well established that the identification of the size-frequency distribution (SFD) of blocks/boulders and the corresponding curve fitting is a pivotal technique that can provide important hints to the formative and/or degradation processes of such features, as showed on the Moon [19][20][21][22], Mars [23][24][25][26][27], as well as on asteroids [28][29][30][31][32] and comets [33][34][35][36]. For this reason, we decided to expand the Martens et al. [1] work analyzing the size-distance relation between the blocks and the fractured Sulci, and obtaining the SFD of the identifiable blocks in the TS area (Figure 1), hence suggesting their generation processes.
The paper is structured as follows: after identifying and counting the blocks located in the Damascus-Baghdad-Cairo Sulci area, we study their SFD, applying a statistical method to evaluate if the data can be fitted by power-law curves. We then discuss the reliability of this approach and study in more detail the block SFD of the three TS, separately. We then discuss our results in the light of similar studies accomplished on other icy studied bodies.

Dataset and Methodology
In order to identify the icy blocks located in the Enceladus SPT, we downloaded the full Cassini ISS-NAC [5] imagery dataset presented in Table 1 of Martens et al. [1], and publicly available through the Planetary Data System (PDS) Archive (https://pds-rings.seti. org/saturn/, accessed on 20 June 2020). Out of the twenty images, we decided to focus on seven (Table 1). This choice was dictated by the fact that such images (i) are not dominated by large shadowed areas that commonly cover more than half of the observed surface, (ii) present a comparable spatial scale, hence favoring the identification of similar block sizes and (iii) cover the widest, almost contiguous area located inside the SPT (Figure 2). This aspect is of great importance since it allows to identify if there are any block SFD differences among the most active TS terrains. In addition to this dataset, we added the ISS-NAC image N1604167158 (indicated in italic in Table 1) that is not present in Martens et al. [1]. This was done to complement the N1597182533 coverage of Damascus Sulcus to the easternmost side. Despite the slightly lower resolution (31.9 m for N1604167158 vs. 25.3 m on N1597182533), this choice allowed a wider block identification, favored by similar illumination conditions. As detailed in Martens et al. [1], we used the Integrated Software for Imagers and Spectrometers (ISIS) in order to convert the raw images into ISIS-readable data cubes, hence attaching all the required Navigation and Ancillary Information Facility (NAIF) SPICE kernels [37] information (mutual position of the target and the spacecraft, camera pointing, as well as spacecraft timing). This was accomplished to map all images into a polar stereographic projection, ready to be imported into the Environmental Systems Research Institute ArcGIS 10.5 software (hereafter ArcGIS).
Through ArcGIS, we then identified the blocks that are located on, or in close proximity to, the TS thanks to the presence of a nearby elongated shadow ( Figure 3A). As commonly done by similar analyses performed on other bodies of the Solar System [26,38,39], once such features have been manually identified, we measured their position on the surface, assumed their shapes to be circumcircles (the maximum resolution of the images does not allow to outline them as polygons) and then calculated their diameter [40], Figure 3B. Finally, to derive the cumulative blocks SFD per km 2 , we computed the corresponding area of the studied surface.
Afterwards, it is possible to evaluate if the resulting SFD is fitted by a power-law or not [41,42]. To do that, we made use of the Clauset et al. [43] method to validate the existence of the power-law fitting model, which is characterized by the scaling parameter called α Cl . We highlight that the power-law index of the cumulative distribution we herafter refer to as "α", is related to the scaling parameter through the following equation: α = 1 − α Cl . The Clauset et al. [43] method also allows the identification of the completeness limit x min , which is the threshold value above which the power-law exists. The estimation of x min is done through the Kolmogorov-Smirnoff (KS) statistic and allows to find the value minimizing it. Afterwards, the parameter α Cl is determined through the maximum likelihood estimator (MLE). The uncertainty for both α Cl and x min is then derived through a non-parametric bootstrap procedure that generates a large number of synthetic datasets from a power-law random generator and performs a number of KS tests to verify if the generated and observed data come from the same distribution. This technique returns a p-value that can be used to quantify the plausibility of the hypothesis. Considering the significance level of 0.10 [43], if the p-value is ≥0.1, then it is possible to conclude that any difference between the empirical data and the model can be explained with statistical fluctuations. On the contrary, if the p-value is <0.1, then the data set does not come from a power-law distribution, but instead, from a different one.

Results
Over the full study area (yellow outline in Figure 2) we identified 17,070 blocks ( Figure 4). The resulting histogram of all counts is presented in Figure 5A. The diameter range of the counted blocks spans between~25-30 m and 366 m, with the biggest block identified~4000 m far from the Damascus Sulcus. As already presented in multiple works, e.g., [22,32], block sizes are only reasonably reliable for sizes larger than 3-4 pixels. Given that the lowest spatial scale of the full dataset is~31 m, only the right tail of the distribution (with diameters larger than 125 m) can be considered meaningful to compute the main statistical properties. The mode of the resulting frequency distribution for diameters larger than 125 m is 133.8 m, the median is 145.3 m with a median absolute deviation (mad) of 21.7 m, while the mean value is 156.4 m, with a standard deviation of 34.9 m. The total considered area is 3296.1 km 2 wide: we can therefore plot the unbinned cumulative number of blocks per km 2 ( Figure 5B) deriving values that span from 4·10 −3 /km 2 at a diameter of 300 m, up to 5·10 −1 /km 2 for sizes~125 m.    By using the above mentioned Clauset et al. [43] methodology, we decided to test if the cumulative distribution function (CDF) can be fitted by a power-law or not. Through the KS statistic we identified that the completeness limit x min for all counts is 149.8 m (with a total number of 688 blocks larger than x min ), while the value of the power-law index α is −5.4 ( Figure 6A). To evaluate the uncertainty for x min and α we generated 5000 synthetic datasets using the non-parametric bootstrap procedure. The scatter-plot of x min against α is presented in Figure 6B in order to evaluate how a change in x min results in a different value for α, while the frequency histograms of x min and α are indicated in Figure 6C,D. The resulting error of x min is 16.2 m, while the one for α is 0.37. As indicated in Figure 6A, the p-value obtained is 0.185. Since this value is ≥0.1, it is possible to affirm that (i) any difference between the empirical data and the model can be explained with statistical fluctuations, (ii) we cannot reject the power-law model for these data, and (iii) the data are consistent with a power-law.   In order to evaluate if there are any possible differences in cumulative size frequency distribution per km 2 or power-law indices, we decided to split the full dataset of blocks into three main groups ( Figure 7): (i) one belonging to the Damascus Sulcus, (ii) one characterized by the blocks situated in close proximity to the Baghdad Sulcus, while (iii) the last one surrounding the Cairo Sulcus. Even if the Baghdad Sulcus branches out in two parts on the studied area, we decided to consider the corresponding blocks as belonging to the same group since the tiger stripe has the same origin.
As done for the full dataset ( Figure 5A), we plotted the histograms of the Damascus, Baghdad and Cairo Sulci blocks ( Figure 8A,C,E), while we plotted the corresponding cumulative number of blocks per km 2 in Figure 8B,D,F. The Damascus area (1208 km 2 ) has a total number of 1838 counted blocks (the lowest spatial scale of the images used for the Damascus block identification is 31.9 m/pixel) and it is the one that shows the largest diameter range, between~40 and 366 m ( Figure 8A). The mode of the skewed distribution for all blocks ≥127 m (4 pixels) lies at 132.4 m, the median at 150.7 m (with a mad of 23.5 m), the mean is 162.0 m, and the standard deviation is 38.2 m. The cumulative number of blocks per km 2 ranges from a minimum of 8.3·10 −4 /km 2 at a size of 366 m to 5.1·10 −1 /km 2 at 127 m ( Figure 8B).
Among the three considered areas, the Baghdad one is the widest with a total surface of 1685. is presented in Figure 6B in order to evaluate how a change in xmin results in a different value for , while the frequency histograms of xmin and are indicated in Figure 6C,D. The resulting error of xmin is 16.2 m, while the one for is 0.37. As indicated in Figure 6A, the p-value obtained is 0.185. Since this value is ≥0.1, it is possible to affirm that (i) any difference between the empirical data and the model can be explained with statistical fluctuations, (ii) we cannot reject the power-law model for these data, and (iii) the data are consistent with a power-law. In order to evaluate if there are any possible differences in cumulative size frequency distribution per km 2 or power-law indices, we decided to split the full dataset of blocks into three main groups (Figure 7): (i) one belonging to the Damascus Sulcus, (ii) one characterized by the blocks situated in close proximity to the Baghdad Sulcus, while (iii) the last one surrounding the Cairo Sulcus. Even if the Baghdad Sulcus branches out in two parts on the studied area, we decided to consider the corresponding blocks as belonging to the same group since the tiger stripe has the same origin. As done for the full dataset ( Figure 5A), we plotted the histograms of the Damascus, Baghdad and Cairo Sulci blocks ( Figure 8A,C,E), while we plotted the corresponding cumulative number of blocks per km 2 in Figure 8B,D,F. The Damascus area (1208 km 2 ) has a total number of 1838 counted blocks (the lowest spatial scale of the images used for the Damascus block identification is 31.9 m/pixel) and it is the one that shows the largest diameter range, between ~40 and 366 m ( Figure 8A). The mode of the skewed distribution for all blocks ≥127 m (4 pixels) lies at 132.4 m, the median at 150.7 m (with a mad of 23.5 m), the mean is 162.0 m, and the standard deviation is 38.2 m. The cumulative number of blocks per km 2 ranges from a minimum of 8.3•10 −4 /km 2 at a size of 366 m to 5.1•10 −1 /km 2 at 127 m ( Figure 8B).
Among the three considered areas, the Baghdad one is the widest with a total surface of 1685.1 km 2 . Out of the 11,883 blocks counted (the lowest spatial scale of the images used for the Baghdad block identification is 25.3 m/pixel), the largest one has a diameter of 343 In order to test if the cumulative distribution functions can be fitted by a power-law or not, we decided to apply the Clauset et al. [43] methodology also to the Damascus, Baghdad and Cairo separate counts. As done for the total counts' analysis, we then evaluated the uncertainty for x min and α generating 5000 synthetic datasets using the non-parametric bootstrap procedure. The results are presented in Figures 9-11, outlined with the same fashion of Figure 6.
For the Damascus case we obtained a x min value of 149.6 ± 14.2 m and a power-law index of −5.11 ± 0.43. The resulting p-value is 0.21. The Clauset et al. [43] methodology applied on the Baghdad counts, instead, returned a x min value of 117.3 ± 16.5 m, a powerlaw index of −5.13 ± 0.34 and a p-value of 0.14. Finally, for the Cairo case we obtained a x min value of 87.8 ± 5.1 m, a α value of −6.32 ± 0.34, and a p-value of 0.11. m and it is located 1150 m far from the Baghdad left branch. The mode of the distribution of blocks ≥101 m is 102.3 m (Figure 8C), the median is 118.5 m with a mad of 19 m, while the mean lies at 127.5 m, with a standard deviation of 30.1 m. The cumulative number of blocks per km 2 at the maximum identified size (~340 m) is 6.0•10 −4 /km 2 and it is comparable to the Damascus one, while at 101 m size the cumulative number of blocks per km 2 is 1.30/km 2 ( Figure 8D).  2.9/km 2 at 74 m ( Figure 8F).
In order to test if the cumulative distribution functions can be fitted by a power-law or not, we decided to apply the Clauset et al. [43] methodology also to the Damascus, Baghdad and Cairo separate counts. As done for the total counts' analysis, we then evaluated the uncertainty for xmin and generating 5000 synthetic datasets using the nonparametric bootstrap procedure. The results are presented in Figures 9-11, outlined with the same fashion of Figure 6.  We highlight that for all three cases the p-values obtained are ≥0.1, hence we can affirm that the use of the power-law curve fitting is meaningful for Damascus, Baghdad and Cairo counts, and any difference between the empirical data and the models are purely explainable with statistical fluctuations. Universe 2021, 7, 82 11 of 18   For the Damascus case we obtained a xmin value of 149.6 ± 14.2 m and a power-law index of −5.11 ± 0.43. The resulting p-value is 0.21. The Clauset et al. [43] methodology applied on the Baghdad counts, instead, returned a xmin value of 117.3 ± 16.5 m, a powerlaw index of −5.13 ± 0.34 and a p-value of 0.14. Finally, for the Cairo case we obtained a xmin value of 87.8 ± 5.1 m, a value of −6.32 ± 0.34, and a p-value of 0.11. We highlight that for all three cases the p-values obtained are ≥0.1, hence we can affirm that the use of the power-law curve fitting is meaningful for Damascus, Baghdad and Cairo counts, and any difference between the empirical data and the models are purely explainable with statistical fluctuations.

Discussion
Given the extremely complex geological setting of the SPT [7] and the ongoing multiple processes (cryovolcanic, tectonic, tidal and thermal) that are occurring in the TS area [6][7][8][9]11,17,18], no single formative process can be invoked to explain the formation, degradation and replenishment of the blocks observed on Enceladus south pole. Nevertheless, before describing which are the blocks' originating processes that can be supported by our analysis, it is worth mentioning which are the ones that can be ruled out. As Martens et al. [1] extensively highlighted, impact cratering cannot be the considered as the driving mechanism originating such blocks, since almost no crater has been found in the SPT.
Seismic disturbance or mass wasting events cannot be invoked either, because they would foresee disproportionately high concentrations of blocks in topographic lows, but the observed icy features appear on top of ridges, slopes and inside valleys (Figures 4 and 7).
On the contrary, by taking into account the vapor density, the drag force, as well as the gravitational acceleration at the surface of Enceladus (0.11 m/s 2 ), Martens et al. [1] evidenced that the ejection process occurring during cryovolcanic eruptions from the TS could be a supported block formative mechanism. It is difficult to conceive that wide enough cracks could let~100 m size blocks pass, since to form such large nozzles both high stresses and strain rates would be needed [44]. However, Hansen et al. [45] indicated that Enceladus collimated jets could reach speeds of 1 km/s or higher, and assuming plume temperatures higher than 296 K, blocks with sizes~100 m could be lifted and ejected [1].
To test the cryovolcanic eruptions as a possible block formation mechanism, for each of our identified blocks we measured the corresponding distance from the nearby TS (Figure 7). In this way we can therefore investigate if there is any block diameter versus TS distance trend, and whether the largest blocks appear closer to the ejection fissures. As presented in Figure 12, when the blocks with sizes smaller than the four pixels identification limit are excluded, it becomes evident that as distance from the TS fissure increases, the blocks' sizes get smaller. Indeed, the smallest blocks (75 m to 130 m size) are present at different TS distance, from tens of m to maximum distances of 20-25 km. On the contrary, the largest blocks (100 s m wide) are generally found at much closer range. This result quantitively supports the ejection hypothesis previously advanced by Martens et al. [1], i.e., where larger blocks usually land closer to the ejection point.
Since the blocks' formation cannot be attributed to a single formative process, another mechanism that needs to be considered is the sublimation of ice in the SPT. Indeed, by modelling the ice sublimation in the 170-240 K temperature range (obtained in the TS area through CIRS data, [46]), Goguen et al. [47] evidenced that sublimation can erode m-scale layers of ice within weeks. Such process is particularly strong close to the geyser vents, i.e., where thermal emissions are concentrated, thermal gradients are larger and hence sublimation rates are higher, but it can still act on the full SPT over longer periods of time, disaggregating blocks from exposed, fractured ice. For this reason, the sublimation process can largely contribute to the evolution of blocks located on SPT. Such phenomenon has not been solely hypothesized on Enceladus [1], but also on another icy body, yet characterized by different size (4 km versus 500 km) and sublimation source (rapidly changing solar insolation versus tidal stress), such as comet 67P/Churyumov-Gerasimenko, herafter 67P.
Thanks to the high-resolution images acquired by the OSIRIS camera [48] onboard the ESA/Rosetta spacecraft, comet 67P is the only case of an icy surface where several boulders/blocks' fields located on different morphological settings have been identified and analyzed. For instance, Pajola et al. [33] studied the cometary protruding blocks located in the heavily fractured, highly sublimating [49,50] Hathor cliff, which is a 900-m high, 1.8 km 2 wide scarp located in 67P small lobe [51]. Besides gravity and tidal stresses, this location can be considered somehow geologically similar to the Enceladus SPT due to its ubiquitous fracture pattern and the active jets, hence allowing the comparison between the two bodies. By studying the Hathor non-dislodged blocks SFD in the 7-40 size range, Pajola et al. [33] derived a power-law index of −5.2 + 0.5/−0.6 (see Figures 11 and 12 of Pajola et al. [33]), which has been attributed to both sublimation and fracturing processes. This value is similar to the one we obtained for the SPT full counts ( Pajola et al. [33] derived a power-law index of −5.2 + 0.5/−0.6 (see Figures 11 and 12 of Pajola et al. [33]), which has been attributed to both sublimation and fracturing processes. This value is similar to the one we obtained for the SPT full counts (  Instead, the Cairo blocks SFD is characterized by a steeper power-law index ( = −6.3 ± 0.3) when compared to the Damascus and Baghdad case. As evidenced in Schröder et al. [32], even for several hundreds of boulders, a difference of 1 in the power-law index may arise simply by chance. This means that the Cairo SFD could still be comparable to the Damascus and Baghdad ones, in particular because the error bars are not small enough to argue this forcefully. Nevertheless, if the Cairo SFD is different from the other two, such steep power-law index could be indicative of a high degree of fragmentation [33]. A similar power-law index has been obtained on locations different from Hathor on comet 67P, Figure 12. Diameter versus TS distance plots obtained for all counts (A), for the Damascus blocks (B), for the Baghdad counts (C) and for the Cairo blocks (D). The vertical green dashed lines are the four pixels limit derived for each case, see Figures 5 and 8. On all plots, the counts located inside the grey shadowed area are not considered because they are below the four pixels identification limit.
Instead, the Cairo blocks SFD is characterized by a steeper power-law index (α = −6.3 ± 0.3) when compared to the Damascus and Baghdad case. As evidenced in Schröder et al. [32], even for several hundreds of boulders, a difference of 1 in the power-law index may arise simply by chance. This means that the Cairo SFD could still be comparable to the Damascus and Baghdad ones, in particular because the error bars are not small enough to argue this forcefully. Nevertheless, if the Cairo SFD is different from the other two, such steep power-law index could be indicative of a high degree of fragmentation [33]. A similar power-law index has been obtained on locations different from Hathor on comet 67P, and it is related to ceiling breakdown, rapid escape of high-pressure volatiles and consequent pit formation, with a resulting blocks' field at its bottom [52,53]. Through the current imaging resolution, we have not identified an unusually higher density of pits inside Cairo, nor it appears to show more erupting vents, or jets (Figure 2, [10]), therefore we cannot invoke such formative mechanism for Cairo's blocks. Such steep SFD may be alternatively related to a locally higher degree of tectonic disruption of lithospheric ice [1] that could be favored by thermal anomalies and ice crust thickness differences that appear to be prominent both inside the Sulci area, as well as at 30-50 km distance from the active SPT location [54]. The large localized thermal gradients may then result in a larger presence of smaller size blocks (hence affecting the steepness of the SFD), when compared to the largest ones. Nevertheless, since Cairo's morphological setting shows no specific peculiarities with respect to Damascus or Baghdad TS (with the present resolution imagery dataset), this point remains to be justified.
Our results suggest that the comparison between Enceladus TS area and Hathor on 67P is intriguing, since it supports the hypothesis that at different size scales, the icy blocks observed on the two bodies could be the result of the coupled fracturing and sublimation activity. This analysis has been made possible because for both cases we have high enough resolution images to identify and measure the sizes of blocks/boulders on their surfaces, hence allowing to advance hypotheses on their formation.

Conclusions
Through the Cassini/ISS-NAC high-resolution images of the Enceladus Damascus-Baghdad-Cairo TS area, we have identified 17,070 blocks with sizes ranging from~25-30 m to 366 m. Then, we identified the blocks SFD and its curve fitting to provide new insights into the formation and/or degradation processes characterizing these features. By using the Clauset et al. [43] methodology we derived that the TS blocks SFD can be fitted by a power-law curve, with a power-law α of −5.4 ± 0.4. The same methodology applied on the three TS separately studied returned a power-law index of −5.1 ± 0.4 for the Damascus case, −5.1 ± 0.3 for Baghdad TS, while a power-law index of −6.3 ± 0.3 for the Cairo area. We explained these results considering that different processes concur in the formation and evolution of such blocks, in particular sublimation and cryovolcanic ejection mechanisms, as previously hypothesized by Martens et al. [1].
To support our findings, we compared our TS SFDs with comet 67P, the only other case of an icy surface where several boulders' fields located on different morphological settings have been identified so far. In particular, the 67P Hathor cliff, where a ubiquitous fracture pattern with several departing jets have been observed, presents a similar SFD (−5.2 + 0.5/−0.6) to the Damascus and Baghdad ones. The sublimation and fracturing processes explaining the Hathor SFD can be applied also to the TS blocks' SFD, suggesting that such mechanisms may favor similar fractures development in the icy matrix, resulting in comparable block disaggregation. For the Cairo case, the steeper power-law index may be indicative of a high degree of fragmentation that could be instead related to a locally higher degree of tectonic disruption of lithospheric ice favored by thermal anomalies and ice crust thickness differences [1], resulting in a larger presence of smaller size blocks. Nevertheless, since Cairo's morphological setting shows no specific peculiarities with respect to Damascus or Baghdad TS, this different power-law index requires further investigation.
Afterwards, we also investigated the block diameter versus TS distance trend in order to test if the cryovolcanic ejection mechanism suggested by Martens et al. [1] can be a possible block formation process. From our counts, we identified that the smallest blocks (75 m to 130 m size) are present from tens of m to maximum distances of 20-25 km from the TS fissures, while the largest blocks (100s m wide) are generally found at much closer range. This result positively supports the vent eruption hypothesis, which foresees that larger blocks usually land closer to the ejection point.
Our findings support that sublimation and cryovolcanic eruptions could explain the formation and evolution of blocks located in the Enceladus SPT, but further higherresolution data are required to test if such mechanisms could be invoked also for smaller sizes. Indeed, m-or cm-scale images are pivotal to largely advance our knowledge of the icy block/boulder SFD formation, as already done in the last decades [20,24,29,32] on different rocky bodies.
In this context, future space missions dedicated to icy worlds will be fundamental; in particular, the images that will be acquired by the JANUS camera [55] onboard the ESA/JUICE mission [56] and the EIS camera [57] of the NASA/Europa Clipper mission [58] will allow to compare and expand the presented study also to the icy surfaces of Europa, Ganymede and Callisto.