Rheological and Mechanical Layering of the Crust Underneath Thumbprint Terrains in Arcadia Planitia, Mars

In the area of Arcadia Planitia in the Northern hemisphere of Mars, mounds indicating fluid and sediment emissions have been already recognized. Here, we show that through fractal and fracture‐spacing analyses of a large vent population it is possible to infer the mechanical layering of the underlying subsurface. Our work includes the mapping of an entire population of 9,028 vents over an area of 122,000 km2. The analysis of mound distribution at the surface led to the formulation of inferences about the subsurface feeding conduits, and to the identification of three mechanical discontinuities at c. 4–5, c. 14–23, and c. 50–55 km. This evidence matches the mechanical stratigraphy recorded by the InSight NASA mission, and is in agreement with independent previous subsurface global modeling, supporting our conclusions.

at several Martian locations (Ivanov et al., 2014;Oehler & Allen, 2010;Orgel et al., 2015;Skinner & Mazzini, 2009). In the specific case of Arcadia Planitia, TPT morphology has been attributed to Amazonian fluid and sediment resurgences (De Toffoli et al., 2019). The pitted mounds that characterize the region in fact display geomorphological characteristics consistent with them being eruption centers, typical thermal inertia of fine-grained loose materials, and spatial patterns highlighting their connection with subsurface fracture systems (De Toffoli et al., 2019), which would be absent whether the mounds would be the result of as ice-related and surficial processes (De Toffoli et al., 2018).
In this work, we present a complete mapping of the TPT at the study area, highlight their main geomorphological traits and spatial distribution characteristics, and finally investigate the underlying subsurface rheology using an analysis of their rooting fracture networks. These fractures, in fact, could have acted as pathways for fluid percolation during the venting process.

Data and Methods
We surveyed an area of ∼122,000 km 2 and mapped more than 9,000 mounds, including the sample area of ∼12,000 km 2 containing around 2,000 mounds previously analyzed by De Toffoli et al. (2019). To produce a reproducible work, we compiled an extensive map of all mounds and all alignments and other significant morphological constituents of the TPT, and performed cluster and fractal analyses and investigated the alignments and spacings within the TPT. These methods allowed us to test the potential relationships between these surface expressions and the underlying systems of connected fracture networks that could have acted as pathways for fluid percolation (Mazzarini & Isola, 2010).
All the observations were performed on a mosaic of 46 Mars Reconnaissance Orbiter (MRO) Context Camera (CTX) images (6 m/pixel resolution; Malin et al., 2007) produced by USGS Astrogeology service's Pilot data portal (Akins et al., 2009), and additional context information was retrieved from THEMIS thermal inertia (Night-time Infrared images: Christensen et al., 2003) data and the Mars Orbiter Laser Altimeter (MOLA; ∼460 m/pixel; Smith et al., 2001;Zuber et al., 1992) topographical elevation model. Mapping and image analysis were performed with ArcGIS® using sinusoidal projections centered on the study area in order to minimize spatial distortion.

Cluster and Fractal Analyses
The fractal analysis adopted here is a statistical procedure based on the spatial distribution of point features and provides information about the connection between emission points at the surface and the subsurface sources of the ejected materials at depth (Mazzarini & Isola, 2010).
The investigation of the subsurface was performed by the application of cluster and fractal analyses to the entirety of the mapped mound population. We mapped 9,028 individual mounds in Arcadia Planitia, over an area of about 122,000 km 2 . For consistency with previous work (De Toffoli et al., 2018, each feature was distinguished from its neighbors through manual identification , and each mapped point was assigned to the mounds' summit pits. We individually mapped mounds above 100 m in diameter (i.e., characterizing traits visible at CTX resolution) and mapped alignments of single or coalescent mounds as lines. While location of each mound was defined by mapping the central pit, lines were drawn connecting the apical pits of coalescent or adjacent mounds. Where chains of coalescent mounds occurred, the summit pit location identification was supported by the observation of perimetral moats or swellings and circumferential troughs associated to the mounds in plan view.
Fractal analysis investigates the spatial properties of vents or fractures and how they fill the space. The investigation of vents, as in the framework of this study, allows us to collect information regarding the underlying fracture systems since vents represent the position of fractures that reached the surface (Mazzarini & Isola, 2010). We performed the spatial distribution analyses according to Mazzarini and Isola (2010) workflow of clustering and fractal investigation. An additional step was taken to further analyze groups of vents where the landform preservation appeared optimal (i.e., minimum erosion, dust cover, obliteration, etc.). To achieve this, independently from the first clustering technique, we also performed a second clustering methodology and subdivided part of the population into additional clusters by selecting groups of mounds where, from a visual analysis of CTX images, minimal erosion or obliteration had occurred (cf. Figure S1 in Supporting Information S1). Subsequently, self-similar clustering (i.e., fractal) analysis was performed on each cluster. This methodology provides insights into the possible presence of hydraulic connections between mounds at the surface and potential fluid sources at depth via an underlying connected fracture network (De Toffoli et al., 2018;Mazzarini, 2007;Mazzarini & Isola, 2010). When the fracture network is a percolating system, its spatial distribution is fractal (self-similar) and is defined in a specific length-scale bounded by a lower and an upper cutoff (Lco and Uco; Mazzarini & Isola, 2010).
Fractal behavior was analyzed by applying the two-point correlation sum method (Bonnet et al., 2001 and references therein). The correlation sum (C(l)) is defined as C(l) = 2 N(l)/N(N − 1) where N(l) are the pairs of points (i.e., the locations of the mapped central pits) whose distance is less than l. If the fractal distribution exists this correlation is valid: C(l) ∼ l D , where D is the fractal dimension. Self-similarity (i.e., fractal behavior) occurs when a plateau is visible in a Δlog(C(l)/Δlog(l)) versus log(l) plot (Clauset et al., 2009;Mazzarini & Isola, 2010). The edges of the plateau define the Lco and Uco, where the Uco represents the distance between the surface and the depth of the fluid reservoirs (Mazzarini & Isola, 2010; cf. Text S1 and Figure S2 in Supporting Information S1). We adopted error minimization procedures to obtain the best estimate for Uco values. We firstly identified the Uco depth range corresponding to the highest values of R 2 (i.e., correlation coefficient of C(l) ∼ l D ), where R 2 is defined by the least squares method. The error of the measure was estimated by calculating half the difference between the maximum and the minimum depths (i.e., distance from the mapped surface) where the R 2 shows the best-fit value (De Toffoli et al., 2018). Additionally, for the estimates to be statistically meaningful, clusters of at least 50 samples are required to extract robust parameter estimates (Clauset et al., 2009). For large datasets (>200) removing a random sample of 20% of the vents does not significantly affect the estimate of the fractal exponent (less than 0.01% of variation) and the error introduced in the estimation of the cut-offs is less than 1%-2% (Mazzarini & Isola, 2010). It has been shown that adding random errors up to 500 m to points locations does not change the estimation of the cut-offs, as they are not statistically different from those computed for the original data set (Mazzarini et al., 2013). In our data set, the number of samples per analyzed cluster span from a minimum of 454 to a maximum of 2,073.

Fracture Spacing Analysis
Fracture spacing analysis is a statistical procedure which allows the characterization of the fracture spacing in a fracture set, defined as a set of fractures of similar orientation and genetic origin. Fracture spacing can be considered as a proxy for the thickness of a brittle layer fractured under tension .
We considered the mound alignments as expressions of fractures that showed tensional behavior, at least during the opening of the conduits that we suggest have allowed the fluid expulsion. We mapped 984 line features over the study area. Each line was mapped by connecting the apical pits of coalescent or adjacent mounds sitting along the same alignment. The phenomenon that this technique investigates is linked to the brittle response of rocks undergoing tensional stress. When layered materials experience extensional stress intense enough to begin cracking, joints perpendicular to the direction of the extension open . The fracturing process evolves from an initial state, when zero joints are present in the system, to a state when it reaches so called "fracture saturation" which occurs when no new joints form and any additional strain is accommodated by further opening of the existing joints . At saturation, joints are organized according to a regular spatial distribution (approximately constant spacing between fractures) and the parametric distribution of spacing is close to a symmetrical normal distribution. In contrast, when fractures are randomly distributed, a parametric exponential distribution of spacing will be observed (Bistacchi et al., 2020). Intermediate situations would show asymmetrical gamma or log-normal parametric distributions (Rives et al., 1992) representing a certain degree of fracture organization evolving toward a situation of fracture saturation Bistacchi et al., 2020;Rives et al., 1992;Tan et al., 2014). Overall, an evolution from strongly skewed to symmetric distributions is observed with the increasing maturity (or saturation) of a joint set.  observed that at saturation, joint spacing tends to be of the same order of magnitude of the thickness of the fractured layer (fracture spacing/layer thickness ratio 1). Thus, the spacing distribution of joints can be considered a proxy for (a) the degree of evolution of the system (shape of the statistical distribution) and (b) the thickness of the fractured layer. To collect spacing data, we did not use all the mapped line features, but selected portions of the study area where the totality of the lineaments was visible at the mapping resolution. This was necessary to avoid errors induced by the obliteration of fractures (e.g., obliteration due to impact debris overlay, dust deposition, lava flows and erosion), that tend to increase the apparent spacing (e.g., in a sequence of three fractures, if the middle one is hidden, the resulting spacing will be twice the real one). We selected eight regions (O1-O9. O2 resulted too small to be statistically significant and it was not possible to run any meaningful analyses. It was accordingly discarded, cf. Data Set S2 in Supporting Information S1) that fulfilled this condition of complete representation.
Since spacing is defined as the distance measured perpendicular to fractures/joints, and since TPTs are sinuous in map view, we automatically computed sets of sampling lines perpendicularly cross-cutting the TPT linear features within each region (code, figures, and complete data set Data Sets S2 and S3 in Supporting Information S1) and measured the spacing between their intersections with the TPT linear features. The sampling lines were obtained as "streamlines" perpendicular to the TPT linear features at every location. The seed points for individual streamlines were manually selected in order to avoid sampling lines with no geomorphological or statistical significance (e.g., extremely short sampling lines or large groups of sampling lines sharing common segments). Preliminary statistics were calculated for single sampling lines in order to test whether the analyses were well posed (e.g., if the region was large enough to yield significant results). When a minimum of 10 exploitable sampling lines for each region was found, we calculated cumulative statistics to combine the sampling lines into a single one, in order to obtain very long sampling lines with a large number of intersections on which to base the fracture spacing. At this stage, empiric spacing distributions can be measured for each individual region, and the parameters of best-fit parametric distributions can be calculated and used as proxies as explained above (Bistacchi et al., 2020).

Observations
TPTs appear as arcuate curvilinear alignments of pitted mounds with a limited relief (visible at MOLA resolution), a summit pit, a lighter relative brightness compared to the surrounding plain, perimeter moats and swellings, a rough surface, and concentric troughs surrounding the central pit ( Figure 2). On Nighttime Infrared (IR) THEMIS imagery data TPTs appeared easily recognizable as dark linear and arcuate structures representing a common low thermal inertia (Farrand et al., 2005;Mellon et al., 2000;Oehler & Allen, 2010). In this area, such mounds have been interpreted as emission centers produced by subsurface sediment mobilization, venting on a large portion of the area under study (De Toffoli et al., 2019;Farrand et al., 2005;Oehler & Allen, 2010;Skinner & Mazzini, 2009;Skinner & Tanaka, 2007;Tanaka et al., 2003). TPTs thus represent highly fractured zones characterized by a dense presence of venting spots. The new mound population herein investigated includes the previously analyzed population (De Toffoli et al., 2019) and mounds in adjacent neighbourhoods, generating a total of 9,028 mapped mounds.
The TPTs appear as arcuate lineaments with a maximum length of 53 km. The alignments are also highlighted in places by surface bulging (Figure 1c). These bulges are observed more frequently in the northern part of the area and almost disappear toward the south. In the study area the bending direction of TPT alignments appear as a prominent trait of the population. In fact, the convexity of lineaments gently bends over the entire population from a E-W toward a N-S direction moving southward (Figure 1a). TPT arcuate alignments show a peculiar curvilinear trend running upslope and indicate interactions with the underlying topography, bending and disappearing where encountering topographical features characterized by higher elevation (Figure 1a).
It has already been demonstrated that a large sample of the mound population is the products of sedimentary volcanism (De Toffoli et al., 2019). At the CTX resolution, we do not observe any perceivable difference between the mounds included and excluded from the sample population, they share homogeneous key traits with very little variability. Accordingly, there is no geomorphological evidence to not assign a homogeneous interpretation to the whole mound field. Overall, their homogeneous morphology and this consistent pattern of behavior leads us to consider the TPT population to be the product of a single process, or a set of processes, that generated non-randomly distributed mound alignments over a defined area.

Rheological and Mechanical Layering of the Crust
The investigation of fractures that underlie the mounds allows a better understanding of processes and conditions leading to the venting of sediments toward the surface as well as assisting the understanding of the subsurface structure. These piping systems connect the surface to the subsurface source of the erupted materials which, in sedimentary volcanism, represent layers enriched in fluids and/or loose sediments that thus constitute mechanical discontinuities. To probe the structures beneath the TPT, we hence applied fractal and fracture spacing analyses that provide information about the rheological boundaries of the underlying rocks (e.g., Mazzarini & Isola, 2010).
Fractal analysis highlighted three main recurrent source depth ranges beneath the TPT field: from the shallowest to the deepest, (a) 4 ± 0.2 km-4.8 ± 0.2 km; (b) 14 ± 2-23 ± 3 km; and (c) 50 ± 5-55 ± 5 km (Figure 3, cf. Table S1 and Data Set S1 in Supporting Information S1). Fracture spacing analysis suggests that the joint network is evolving toward saturation (we observed log-normal and normal distributions) and suggests an uppermost mechanical layer thickness of the same order of 4 km as the shallowest material source detected with fractal analysis (Figure 3). Arcadia Planitia is a sedimentary basin in the Martian lowlands (Tanaka et al., 2014) and the first few kilometers from the surface are likely composed of a sedimentary sequence of basin infilling with a possible partial surficial cementation due to shallow sulphate-bearing ice or brines (e.g., Cooper & Mustard, 2002;Golombek et al., 2020). This sediment pack has the potential to also represent the source of loose materials that were vented out from the mound pits in the TPT.
Deeper, we detected an intermediate source ranging from between 14 and 23 km depth. Noticeably, this interval matches the estimated range of depths of the proposed clathrate-rich cryosphere stability zone up to the proposed cryosphere/hydrosphere contact (Clifford et al., 2010) leading to the inference that this fracture network system had the potential to connect the source of the ejected fluids (ostensibly water resulting from ice melting or gas hydrate dissociation) to the surface. Hence, given the observed layering, we propose that the fluids were probably sourced from the deeper level and mixed with loose materials once they reached the topmost layer, which acted as an intermediate reservoir before the final discharge.
Moving deeper still, fractal analysis suggested the presence of a deep source layer located between 50 and 55 km below the surface. To crack such an extensive portion of the crust, volcanic, tectonic or impact processes must be invoked (Rossi, 2018). A single impact event, as observed in Arabia Terra , can be excluded since there is no evidence on the surface of an impact large enough to generate a 122,000 km 2 mound field. To include the entire field, a crater of over 500 km radius should be present. Such extensive dimension is reached by major basins like Argyre, while in Arcadia there is no evidence at the surface of an impact of this magnitude. However, a combination of multiple impacts, today buried and obliterated, might have had the potential to deeply fracture the crust, generating organized fracturing that  (Clifford et al., 2010). In blue, the brittle-ductile transition for wet rheology under present-day (in orange) and 1 Ga (in red) conditions (Azuma & Katayama, 2017). (c) In the furthest right column, the reconstruction of the rheological stratigraphy inferred from fractal and fracture spacing analyses at the study site is shown.
pre-existed the TPT (e.g., arcuate patterns). Tectonic processes might also be responsible for such pervasive deep cracking. The brittle lithosphere (i.e., schizosphere) of recent Martian times (∼1 Ga) is assumed to reach depths comparable to those herein observed (Azuma & Katayama, 2017;Grott & Breuer, 2008 and could have been fractured by prolonged tectonic activity Lognonné et al., 2019). In this context, volcanic processes might also have played a role. Dark mounds that might be of magmatic origin displaying possible summit calderas (Figures 2d-2g) are also present in the area, but they appear clearly distinct from the mounds included in our analyses. Additionally, in some cases the TPT mounds appear to overlap such features (Figures 2d and 2f). In view of this, we do not favor the hypothesis that different groups of mounds (i.e., TPT) may have originated from different processes. Indeed, the brittle-ductile interface (i.e., rhizosphere-plastosphere boundary) is a rheological discontinuity which can act as a boundary where magmatic, hydrothermal and epithermal fluids can accumulate and eventually reach the surface when favorable cracking occurs.
Overall, the application of fractal and spatial analyses allowed us to identify mechanical discontinuities (i.e., materials reservoirs and brittle layer thickness) and thus reconstruct the mechanical stratification in the study area (Figure 3), providing additional support to the mechanical discontinuities of the Martian subsurface separately reported by various authors (e.g., Azuma & Katayama, 2017;Clifford et al., 2010;Giardini et al., 2020;Lognonné et al., 2020).
Martian mechanical stratification was also inferred by the NASA mission InSight (Interior Exploration using Seismic Investigations, Geodesy and Heat Transport) Lognonné et al., 2020), the first robotic mission aimed at exploring seismicity on planets other than Earth. InSight recorded seismic activity on Mars due to local perturbations, such as wind, diurnal temperature changes, and tectonism with spectral features very similar to those observed for seismic events on Earth Golombek et al., 2020). In addition, for the first time a local mechanical stratigraphy for a sector of the planet has been inferred by analyzing the propagation of the seismic waves Lognonné et al., 2020). At the InSight landing site, the stratigraphy appears to be characterized by a very thin and weak surficial regolith layer (cohesion up to 12 kPa, and Young's modulus of 0.2 GPa; Golombek et al., 2020) 10-100 m thick ; an upper crust layer 5-7 km thick probably composed of regolith, sediments and highly altered and damaged basalts (with shear-wave velocity V s 1.7-2.1 km/s) and a deeper crust layer, from 5 to 7 km to a depth of about 15-20 km, composed mainly of less altered basalts (V s 2-3 km/s; Lognonné et al., 2020). Beneath this layer, down to 100 km depth, the crust is composed mainly of basaltic rocks (V s 2.8-4 km/s) .

Wider Implications
On Earth, the thickness of different mechanical layers (i.e., sedimentary basins and brittle crust) control the scaling laws of fractures and earthquakes (e.g., Davy, 1993;Ouillon et al., 1996;Pacheco et al., 1992). In volcanic areas, the spatial distributions of both faults and vents provide a coherent image of the crustal layering as observed in the areas along the East African Rift System (e.g., Mazzarini & Isola, 2010. The fractal distribution of mounds/vents correlates well with the crust layering and structure providing the depth of the fluid source for mud volcanoes in the Greater Caucasus region in Azerbaijan (Bonini & Mazzarini, 2010). The fractal analysis of spatial distribution of vents has also been successfully applied to other areas of Mars  as well as on icy planetary bodies, such as Enceladus and Ganymede (Lucchetti et al., 2017(Lucchetti et al., , 2021. In view of the homogeneity of characterizing traits and structure distribution patterns of TPTs, the TPT constituting mounds of Arcadia Planitia are likely to have a common origin. Subsurface sediment mobilization is usually the result of the upwelling of fluids stored at depth that occurs when the buoyancy forces that push the mixtures upward exceed the confining lithostatic pressure (Kopf, 2002). This can be triggered by fluid supply, viscosity decrease or increased pore pressures that in turn could be induced by various events such as impact cratering unloading, hydrothermal pulses, tectonic loading, seismic triggers or injection of gas following clathrate dissociation (Dimitrov, 2002;Oehler and Allen, 2010;Prieto-Ballestreros et al., 2006;Skinner & Mazzini, 2009;Skinner & Tanaka, 2007). Since no unequivocal evidence distinguishing these phenomena was found from the analyses of the orbital optical imagery (e.g., faults overlapping, large volcanic/mud flows, etc.), the trigger that led to the generation of the TPTs could not be identified.
In relation to this question, we report the pervasive presence of linear structure that we interpret as faults in the study area, since this kind of feature has the potential to play a role in subsurface mobilization, given their involvement in some of the above-mentioned processes such as tectonics, liquefaction or hydrothermalism. We observed a set of N-S fault scarps at least 100 km long, characterized by ∼60 m of topographic drop, around and intersecting the whole TPT population (with TPTs overlapping in places the fault traces) and located 20-90 km apart from one another (Figure 4b). It needs to be taken into account that no relationship between large-scale faults and liquefaction events has ever been previously observed on Mars. In addition, there is no proof that the observed faults were active while water/ice reservoirs were available in the Martian subsurface. Nevertheless, on Earth this connection is well established, and the link between faults and fluid-sediment resurgence is commonly observed, since material upwellings preferentially develop along pre-existing surface or near-surface discontinuities and weaknesses (e.g., Bonini et al., 2016;Manga & Brodsky, 2006;Manga et al., 2019). Excellent examples of this phenomenon are the liquefaction events related to the 2,010-2,011 Canterbury earthquake sequence in New Zealand (Quigley et al., 2013;Reid et al., 2012;Townsend et al., 2016). In this case, liquefaction occurred in areas (up to ∼70 km from the epicenter) of abandoned or recently in-filled river channels, swampy ground and areas underlain by estuarine deposits; the most heavily affected areas appear to be related to near-surface alluvial features, such as channels and meander loops, filled with loose, unconsolidated sediment (Bastin et al., 2015;Cubrinovski et al., 2011;Townsend et al., 2016;Villamor et al., 2016;Wotherspoon et al., 2012). In this context, sand boils and sand volcanoes developed aligned along the pre-existing buried features, generating arcuate or straight sets of discharge vents alignments (Reid et al., 2012;Townsend et al., 2016).
On Mars, mid-latitude zones between 30° and 55°, including the Arcadia region, recorded traces of glacial features up to the late Amazonian period, such as glacier-like forms and viscous flow of ice-rich sediments capable of moving significant volumes of material, that may indicate the presence of buried moraines at the study site (Brough et al., 2019;Souness & Hubbard, 2012). In this particular case, the TPT convolute geometries might reflect the arrangement of buried moraines (e.g., Arfstrom & Hartmann, 2005) which may have provided the pre-existing points of weakness (i.e., unconsolidated debris) finally exploited as preferential pathways by fluids and sediments to discharge to the surface after their remobilization (Figure 4a). This process might have taken place where fluids and fractures reached the surface as a result of the deeper subsurface mechanisms and rheology, which determined fracture and vent spacing.

Conclusions
The investigation of planetary interiors and subsurfaces is a fast-developing field, tackling largely unknown environments. On Mars, thanks to the InSight mission, a new era for planetary seismology and geology has opened, and many more discoveries are yet to come. We carried out independent analyses on chains of numerous mounds interpreted as sedimentary volcanoes, retrieving the mechanical stratigraphy for Arcadia Planitia, and compared the results obtained with the observations performed at the InSight landing site, and multiple subsurface modeling studies. Through the application of fractal and fracture-spacing analyses, we identified three mechanical discontinuities at c. 4-5, c. 14-23, and c. 50-55 km, which are in good accordance with the layering inferred from the seismic wave propagation recorded by InSight, which revealed an upper weaker layer at 4-7 km over a bedrock up to 15-20 km in depth. Our outcomes also match previous global modeling of the Mars interior that predict a cryosphere base depth at c. 18 km (Clifford et al., 2010) and a brittle lithosphere base depth comparable to the deep discontinuity observed in the present study at c. 50-55 km (Azuma & Katayama, 2017;Grott & Breuer, 2008. Further modeling and observations are needed to tackle and explain the processes that triggered, or contributed to, the formation and developed of the fracture systems and materials upwelling. An interpretation investigating the origin of the mounds might in fact provide further clues on recent water-related processes in the Martian crust.
By demonstrating the efficiency of these techniques in the field of subsurface mechanical layering reconstruction, we provide a workable tool to expand these investigations to further areas on Mars. This has the potential to significantly increase the number of case studies beyond the range of InSight, but that can assist in a global mapping of the subsurface layering of the planet.