Numerical Simulations of Drainage Grooves in Response to Extensional Fracturing: Testing the Phobos Groove Formation Model

A layered heterogeneous structure is adopted for Phobos to predict its tidal response as shown in Figure 1(b). The Phobos interior is represented as a weak core, barely holding together as a mechanical rubble pile, surrounded by layers of regolith represented by discrete particles. This highly damaged, weak interior structure is consistent with density data in Andert et al. (2010) and gravity data in Guo et al. (2021). The weak interior allows stress from the increasing tidal field to readily deform the interior, thereby forcing the outer regolith to reshape and experience surface ruptures (Hurford et al. 2016). Note that the friction forces between constituent gravels can provide enough structural strength for this cohesionless granular body to remain intact inside the fluid Roche limit. The outer regolith consists of a cohesive subsurface layer (100 m) and an upper loose layer (50 m), consistent with structural interpretations of asteroids Ryugu (Arakawa et al. 2020) and Bennu (Barnouin et al. 2019).

Phobos is tidally evolving on a timescale of tens of millions of years, based on the exchange of angular momentum with Mars. The orbital migration dynamics of Phobos is derived in Appendix A, which, combined with its possible initial positions, can be used to reconstruct the orbital history of Phobos since its origin (see Appendix B). For the proposed layered structure for Phobos, we adopt a thin-shell spherical approximation as suggested in Kattenhorn & Hurford (2009) and Hurford et al. (2016) to calculate the surface tidal stresses exerted by its evolving global shape due to orbital decay (see Appendix C). The evolving stress field is then applied to the surface regolith to predict the response of Phobos driven by its tidal deformation.

To explicitly model the particle dynamics of granular regolith during the tidal deformation processes, we perform discrete-element simulations using an original soft-sphere discrete-element model code DEMBody (Cheng et al. 2019, 2020), which is capable of accurately simulating the complex behavior of particle–particle and particle–boundary interactions. In this code, the nonlinear contact forces between two contacting objects (sphere–sphere, sphere–wall) are calculated based on Hertz–Mindlin–Deresiewicz contact theory (Somfai et al. 2005) according to the particle overlaps (typically <1% of their reduced radii). The implementation of SSDEM in DEMBody consists of four components in the normal, tangential, rolling, and twisting directions, which indicate the elastic forces, damping resistances, and frictions in contact point (Cheng et al. 2018). A friction coefficient μ is used to control the stick-slip Coulomb friction between colliding particles, and a shape parameter β, which represents a statistical measure of real particle nonsphericity, controls the rotational resistance to relative motions in rolling and twisting directions (Jiang et al. 2015). A second-order leapfrog integrator with a time step of 0.02 s is used for integration of the equations of motion of the granular system. Note that the shear failure of a granular system is triggered by microscopic sliding or rolling between particles in contact; therefore, the combination of the above two friction coefficients controls the overall friction angle of granular materials. A dry cohesion model, parameterized by an interparticle cohesive strength c _p , is adopted to mimic the weak attractive force within granular regolith (Zhang et al. 2018). The contact between the two particles will break if the external tensile force along the normal direction exceeds the cohesive force; therefore, the interparticle cohesive strength c _p stands for the macroscopic tensile strength of granular materials.

The mechanical properties of Phobos's surface are basically unknown; however, some remote-sensing observations and laboratory experiments can help constrain the parameter space (Ballouz et al. 2019). Radar albedo data in Busch et al. (2007) imply a near-surface porosity of Phobos's regolith of at least 40% ± 10%. Photogrammetric measurements show that the steepness of local slopes on Phobos's surface can reach 30°–40° (Basilevsky et al. 2014). We thus set μ and β to 0.6 and 0.6, respectively, corresponding to sand particles with medium hardness and material friction angle of ∼30° as determined by laboratory tests (Jiang et al. 2015). The free parameter c _p is used to adjust the material tensile strength T, ranging from c _p = 0 to 48 kPa. Note that c _p is the microcohesion between contacting particles, while T is the macroscopic tensile strength of granular materials.

Considering the expansive parameter space and significant computational expense per run, we restrict our simulations to 23 local regolith patches across the surface of Phobos as shown in Figure 1(b), rather than simulating the global outer layer. We set up simulations of regolith beds with length of 800 m, width of 800 m, and height of 150 m initially. The regolith layers have an initial porosity of 40% and are made up of 3,576,846 grains with material density of 2.8 g cm⁻³ and radii between 1.0 and 2.0 m (Andert et al. 2010). The particle size distribution is selected to allow numerically resolvable groove structures to develop in a reasonable runtime (∼10 days), which is a compromise between simulation accuracy and computational efficiency. Once the granular bed settles down under the Phobos gravity (0.57 cm s⁻²), particles with initial height below 100 m are set to be cohesive with corresponding cohesive strength c _p , and the whole system is then allowed to re-equilibrate to remove possible local disturbances. Deformable periodic boundary conditions (DPBCs) are imposed on the simulated patch in the directions perpendicular to the initial local gravity vector to allow the grains to deform freely in an unconfined environment as suggested by Ballouz et al. (2019). DPBCs are implemented by replicating a patch of particles in the horizontal directions. Each replicated patch contains "ghost" particles that match the relative positions of the original particles. The sizes of the simulated patch in the x- and y-directions are adjusted according to the corresponding deformation velocities. In addition, for continuity, when a particle crosses one of the deformable boundaries, its velocity component in this direction is modified to account for the cell deformation. Therefore, depending on whether it crosses over in the direction of increasing coordinates or vice versa, its velocity component in this direction will be adjusted by ±2 v_deform, where v_deform is the deforming velocity of the patch boundary that the particle crosses.

After this setup process, the simulation patch is expanded and/or shrunk gradually in the directions of local principle stresses, mimicking the deformation of the regolith layer as Phobos's global shape evolves owing to orbital migration toward Mars. For simplicity, the gravity in each patch is assumed to remain constant in the normal direction during deformation.

The deforming velocities v_deform of patch boundaries are determined by the corresponding local tidal stresses, specifically the ratio of principle strain rates derived from the tidal stress field given by the thin-shell model (see details in Appendix A). Figure 1(c) gives an example of the biaxially deformable periodic boundaries used for the modeled patch located at 60° N and 0° E.

As for varying the strength, we used nine different values of cohesion, but only for those of the 23 patch locations that had the most differing stress ratios. Considering the computational cost per simulation, for patch locations with similar stress ratios only one simulation was performed, c _p = 36 kPa, allowing for across-the-board comparison, for a total of ∼100 simulations.

We take the patch located at 60° N and 0° E with c _p = 36 kPa as a representation to detail the fracturing and draining dynamics of the granular beds during the tidal deformation process (see the animation embedded in Figure 1). The regolith patch is forced to reshape, as described above, and accumulates failure stresses on external loading. The cohesive layer deforms elastically at first. When the tensile yield strength T is exceeded, interparticle bonds fail and break, exhibiting quasi-brittle behavior. Parallel fractures are resolved, which initiate with nearly linear traces over a uniform spacing, whose appearance is slightly corrugated and oriented perpendicular to the extension direction. This regular organized pattern resembles the transverse crevasses on terrestrial glaciers under tensile stresses (Colgan et al. 2016), but here for much weaker materials in much lower gravity. As these cracks open up, the failure stresses intensify at the front of the fracture tip (Bassis & Jacobs 2013) and eventually overcome the overburden pressure. Therefore, new microcracks nucleate ahead of the tip and propagate downward perpendicular to the principal tensile stresses, forming a set of subvertical opening mode fissures, as shown in Figure 1(d).

Where fracture openings occur beneath unconsolidated surface material, they cause subsidence of the overlying loose sands into the crevice-generated void space (Horstman & Melosh 1989; Ferrill & Morris 2003; Wyrick et al. 2004), forming separate crateriform depressions that lack elevated rims. These pit craters align above the underlying dilational crack, whereby they develop a quasi-linear pattern naturally. As deformation goes on, particles near the pit walls migrate progressively toward the floor sinkholes in response to the expanding fractures. Initially isolated pits thus gradually grow and intersect, forming an elliptical morphology with the long axis along the chain. In the later stage of the drainage process, adjacent pits further merge into scalloped-sided troughs, and collapsing depressions ultimately coalesce into relatively smooth-sided continuous grooves (Figure 2 shows the development sequence). We analyzed the pattern of model-predicted grooves and found that most of them depart from the plane perpendicular to the tensile direction by less than 50 m over the simulated patch. The geomorphic transitions from pitted depressions to scalloped troughs to linear grooves that are observed on Phobos probably imply different stages of drainage development and are all reproduced in our tidal grooving scenario.

Figure 3 demonstrates the cross-sectional profiles of the grooves formed by granular drainage. Most model-predicted grooves appear shallow. Morphological analysis shows that they are typically 10–30 m deep and 80–120 m wide. Slopes of groove walls vary from 10° to 30°, with an average value of 25°, which are around the angle of repose for loose materials. The simulated grooves also all have flat edges and lack raised rims that are typically associated with impact structures; thus, this could serve as a diagnostic of drainage grooves versus other formation scenarios, such as impact crater chains or boulder-excavated troughs, that will be studied in much further topographic detail using data from future missions. A number of the subsurface fractures show narrow and steeply plunging profiles with vertical or near-vertical walls, morphologically resembling extensional crevasses on terrestrial glaciers (Colgan et al. 2016), but here for much weaker materials in much lower gravity. This suggests that there could be underlying canyons inside Phobos hidden under the observed surface grooves, if they are created by tidal stresses.

Zoom In Zoom Out Reset image size

To consider a range of potential cohesion strengths, we systematically vary the interparticle bonding strength c _p for different simulated patches across Phobos's surface. A summary of granular failure patterns for a range of c _p and patch locations (determining the tidal stress ratio) can be seen in Figures 4 and 5, which reveals well-defined transitions between four regimes of outcomes as distinguished by outline colors: (i) intact deformation for lower cohesion cases, in which the layer distorts without surface fissures, showing high ductility upon mechanical loading; (ii) polygonal failure for higher stress ratio cases (biaxial stretching), in which two or more crosscutting fractures with disparate orientations coexist, sculpting the topography to a complex craquelure-like appearance that is similar to the wrinkling patterns observed on biaxially compressed shells (Al-Rashed et al. 2017); (iii) fragmented forms for intermediate cohesion cases, in which several subtle discontinuous depressions develop, exhibiting scattered segments with lateral offsets; and (iv) parallel patterns such as in Figure 1(d) for higher cohesion cases, showing quasi-brittle behavior during the grooving process. These distinct regimes are strongly dependent on the cohesion strength and the stress state, plausibly due to the plastic-to-brittle transition of cohesive granular materials (Iveson & Page 2004; Smith & Litster 2012; Alarcon et al. 2018). The intriguing patterns reveal the complexity of the failure behavior of granular mediums. The results here could be applied not only to Phobos grooves but also for general understanding on the granular failure mode of small-body regolith under tensile/compressing loading.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Results imply that grooving parallel stretch marks across the midlatitudes of Phobos surface is possible in larger cohesion cases such as c _p = 42 kPa or even larger. This corresponds to a macroscopic material tensile strength of T ∼ 1.0 kPa (see Appendix C for the translation from microscopic parameter to macroscopic quantity of granular matter). We map the prominent linear depressions on Phobos and find that most grooves at midlatitudes line up with the model-predicted groove orientations as shown in Figure 6. This consistency demonstrates the possibility that some observed linear striations might be the surface manifestations of underlying canyons hidden inside the small moon. The low gravity and tenuous cohesion, combined together, make the powdery but mildly cohesive layer show brittle behavior, a characteristic once thought to be reserved for coherent materials. This unique process produces well-defined granular fractures on low-gravity regolith. If true, then many linear features on small bodies (Sullivan et al. 1996; Watters et al. 2011) and icy bodies (Morrison et al. 2009; Spencer et al. 2020) that are previously believed to form through tectonism on monolithic bodies with structural continuity could plausibly be granular landforms with mild cohesion.

Zoom In Zoom Out Reset image size

However, there are significant regions of Phobos showing absences and anomalies compared with the tidal fracturing model predictions. No grooves are found in a ∼50° wide region centered around the trailing apex (Murray & Heggie 2014). In the tidal formation scenario, this might be attributed to a paucity of cohesive material (i.e., coarser regolith; Scheeres et al. 2010), so the tidal strain might create more continuous deformation and not fissures as predicted in Figure 4. Such regional compositional differences have been invoked to solve the weakness stated here by Hurford et al. (2016). This explanation might be consistent with the observation that discontinuous depressions predominate near the absence zone (Thomas et al. 1979; Murray & Heggie 2014), possibly the marks of underdeveloped subsurface fractures. In addition, near the sub-Mars region the tidal fracturing model predicts almost isotropic tensile stresses and polygonal fractures. However, only east–west linear grooves are observed in this location (Murray & Heggie 2014). A plausible explanation is that stress concentrations (Bruck Syal et al. 2016), especially the long-wavelength topography associated with Stickney (Hurford et al. 2016), make the stress field anisotropic compared to a spherical model. Crisscrossed grooves cluster in the equatorial zones of the leading side. One subset is still aligned with the current orbital decay stress field. The other one, however, crisscrosses the tensile stresses. This mismatch was associated with possible multiple stages of tectonic extension and direction by Hurford et al. (2016). It means that ancient stress fields of Phobos that are different from its current stess fields generated the mismatched grooves.

Other mechanisms are also thought to have contributed to groove formation, such as rolling and bouncing boulders (Ramsley & Head 2019) and chains of impact craters formed by ejecta from Mars (Murray & Heggie 2014) or Phobos (Nayak & Asphaug 2016). There is a possibility that more than one mechanism is responsible for groove formation on Phobos, and that the tidal fracturing model is one of them, or the correlation between the tidal stress orientations with Phobos grooves in midlatitudes is just a coincidence. It is premature to draw conclusions beyond the scope of what is presented here. A more detailed study including the 3D shape of Phobos, as well as studying these absences and anomalies of grooves, would be a suitable application of the higher-resolution topographic data anticipated from the Mars Moons eXploration (MMX) mission (Kuramoto et al. 2022), which may conclusively determine the origin of these enigmatic striations.

Assuming that the tidal fracturing model is at work on Phobos, the tidal stresses accumulate as Phobos drifts inward, sculpting the striking grooves on its surface. These linear striations, therefore, serve as "tree rings" recording the past tidal evolution and surface modification on the Martian moon. This makes the orbital history a powerful tool, when combined with the grooving scenario, for inferring the hidden structure and formation history of Phobos. Competing theories suggest either that the Martian satellites are captured asteroids (Pajola et al. 2013) or that they were formed following a giant impact on Mars (Citron et al. 2015; Hesselbrock & Minton 2017; Canup & Salmon 2018). The extrinsic origin, however, is difficult to reconcile with the current nearly circular and coplanar orbits of Phobos and Deimos (Canup & Salmon 2018). Alternatively, this would be a natural consequence of the moons coalescing from an ancient equatorial disk around the red planet. A basin-forming impact on Mars over 4.3 Gyr ago has been invoked for the formation of the debris ring that proceeds to accrete Phobos- and Deimos-like moons directly, i.e., a single giant impact model (Citron et al. 2015; Canup & Salmon 2018). However, such an early forming of Phobos requires an initial orbital height far from Mars (Rosenblatt et al. 2016; Hesselbrock & Minton 2017), near the synchronous orbit. This means that Phobos should pass a 2:1 resonance with Deimos as it moved inward, which would probably excite the eccentricity of the outer satellite to an order of magnitude greater than its current value (Yoder 1982; Ćuk et al. 2020). A cyclic Martian ring-moon hypothesis (Hesselbrock & Minton 2017; Ćuk et al. 2020) is a possible formation mechanism to solve this conundrum. In this model, Phobos, with a much younger age, is the latest generation of moons born in the repeating ring-forming episodes. The most recent cycle, which produced present-day Phobos, would have then occurred near 3.34 Mars radii as deduced from the equipotential figure of Phobos (Hu et al. 2020).

Based on the two possible in situ accreting formation scenarios, the orbital migration history of Phobos can be reconstructed respectively as shown in Appendix A. We adopt four internal models, with low rigidities of μ_i = 10⁵, 10⁶, 10⁷, and 10⁸ Pa, consistent with its low density, which implies a large porosity broadly (Andert et al. 2010; Taylor & Margot 2011). The cohesive layer is assumed to have a rigidity of μ_o = 10⁸ Pa, to be consistent with Hurford et al. (2016) based on an analog to the Lunar regolith layer. Its past orbits, in conjunction with the layered tidal stress model (see Appendix A), are used to deduce the tensile stress field experienced on the de-orbiting Phobos (Figure 7).

Zoom In Zoom Out Reset image size

We consider the maximum stress criterion to estimate the stress required to fracture the granular regolith, which assumes that tensile cracking occurs once the maximum tensile stress exceeds the corresponding material yield strength T (estimated to be at least ∼1 kPa to reproduce parallel grooves; see Appendices B and C). In the single giant impact scenario, Phobos started its inward migration from about 5.76 Mars radius at 4.5 Ga. Assuming that Phobos has experienced the same solar system impactor flux as Mars (Schmedemann et al. 2014), Stickney is then estimated to be 4.0–4.2 Gyr old through crater statistics. This impact event should predate the groove formation given the stratigraphic order that grooves superpose Stickney; therefore, most grooves probably began forming about 1.5–3.8 Ga ago, corresponding to interior rigidities between 10⁶ and 10⁷ Pa (Figure 7(a)). For Phobos with a moon-ring cyclic origin, in which it evolved from near 3.34 Mars radii at 93 Ma, the fracturing of Phobos's surface probably could date back to 60–90 Ma, with interior rigidities restricted within 10⁶–10⁷ Pa (Figure 7(b)). This more recent estimate of the grooves' age seems to be in tension with the crater-age determination mentioned above. We speculate that while the body was accreting, impacts with ring material and sesquinary ejecta would have contributed to the overabundance of craters (Nayak & Asphaug 2016; Hesselbrock & Minton 2017). Such a planetocentric bombardment, mixed up with background cratering, makes it easy to overestimate the crater age.

As a result of tidal dissipation, Phobos, which orbits within the synchronous limit, evolves gradually inward. We made a reasonable assumption of tidal locking in 1:1 spin–orbit resonance during its orbital decay (Black & Mittal 2015; Samuel et al. 2019). By considering the exchange of angular momentum with Mars, the evolution of the satellite's orbit (Charnoz et al. 2010; Hesselbrock & Minton 2017) is given as

$$\begin{eqnarray}&&\frac{da}{dt}=-\frac{3{k}_{2}{{nM}}_{\mathrm{Ph}}{R}_{{\rm{M}}}^{5}}{{{Qa}}^{4}{M}_{{\rm{M}}}}\left[1+\frac{51{e}^{2}}{4}\right].\end{eqnarray} \tag{ A1 }$$

Here k₂ = 0.148 ± 0.017 is the degree-2 tidal potential Love number of Mars and Q = 88 ± 16 is its tidal quality factor (Nimmo & Faul 2013). Parameters n, e, and a are the mean motion, eccentricity, and semimajor axis of Phobos, respectively. M_Ph and M_M are the masses of Phobos and Mars, and R_M is the radius of Mars. Because the orbital eccentricity of Phobos is very small (∼0.0151), its eccentricity change is neglected (Charnoz et al. 2010; Hesselbrock & Minton 2017).

The origin of the Martian satellites, Phobos and Deimos, is not well understood. An external origin has been proposed, in which the two moons are suggested to be captured asteroids (e.g., Pajola et al. 2013); this scenario, however, does not readily explain the current small inclination and eccentricity of Phobos and Deimos (Burns 1978). Other scenarios advocate for an in situ origin from an ancient equatorial disk around the planet, such as might be created by a single giant impact (Citron et al. 2015; Canup & Salmon 2018). In this hypothesis, Phobos formed directly from the coalescence of a debris disk ejected during a violent collision with the Martian surface over 4.3 Gyr ago, for instance, the impact responsible for Hellas basin (Craddock 2011) or the planet's notable hemispheric dichotomy (Marinova et al. 2008). The challenge of disk coalescence models is that the impact-generated disk must be massive enough to be gravitationally unstable and not disperse, making it difficult to form such small satellites.

For parameters of k₂ = 0.148 and Q = 88 for Mars, integrating back in time from Phobos's current orbit at 2.76 R_M implies that at 4.5 Ga the orbital semimajor axis of Phobos would have been ∼5.76 R_M. However, such an orbital migration of Phobos would probably result in a 2:1 resonance with Deimos that might push the outer satellite onto a strongly elliptical orbit, which is inconsistent with its present state (Hesselbrock & Minton 2017). A possible solution to this conundrum, which could also explain the small mass of Phobos and Deimos, is that Phobos was originally a much larger body that coalesced from a massive debris ring produced in a giant impact. Forming inside the corotation radius of Mars, it then experienced cycles of progressive tidal breakup (Hesselbrock & Minton 2017; Ćuk et al. 2020). In this cyclic Martian ring-moon hypothesis, Phobos is the latest in a sequence of in-spiraling moons that were tidally disrupted inside the Roche limit. The largest remnant of the disruption (currently Phobos) would then migrate outward owing to resonant interactions with its own tidally stripped debris. This remnant would then spiral in again owing to tidal dissipation in Mars, repeating the cycle.

Agnostic of the nature of its origin, equipotential analysis of Phobos (Hu et al. 2020) suggests that the current satellite might have been originally emplaced near 3.34 R_M, giving a relatively young age of 93 Myr deduced from its current location. This leads us to consider two end-member candidates for the orbital histories of Phobos, as plotted respectively in Figure 8: on the one hand, beginning at ∼5.76 R_M at 4.5 Gyr, and on the other hand, beginning at 3.34 R_M at 93 Myr.

Zoom In Zoom Out Reset image size

We consider Phobos as having a layered heterogeneous structure, made up of a homogeneous spherical interior and a discrete granular outer layer. Surface stresses are computed using a thin-shell spherical approximation (Kattenhorn & Hurford 2009). The horizontal strain in the elastic shell, induced by the increasingly tidally distorted interior, results in stresses on the surface (Hurford et al. 2016) given by

$$\begin{eqnarray}&&{\sigma }_{\theta \theta }=\displaystyle \frac{9{M}_{{\rm{M}}}{\mu }_{o}{h}_{2}}{88\pi {\rho }_{\mathrm{Ph}}}\left(\displaystyle \frac{1}{{a}_{f}^{3}}-\displaystyle \frac{1}{{a}_{i}^{3}}\right)\left(5+3\cos 2\theta \right),\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&{\sigma }_{\varphi \varphi }=-\displaystyle \frac{9{M}_{{\rm{M}}}{\mu }_{o}{h}_{2}}{88\pi {\rho }_{\mathrm{Ph}}}\left(\displaystyle \frac{1}{{a}_{f}^{3}}-\displaystyle \frac{1}{{a}_{i}^{3}}\right)\left(1-9\cos 2\theta \right),\end{eqnarray} \tag{ A3 }$$

where θ is the colatitude with respect to the tidal distortion's axis of symmetry and σ_{θ θ} is the stress along the surface in the direction radial to the axis of symmetry, while σ_{φ φ} is along the surface orthogonal to σ_{θ θ} (Hurford et al. 2007). Parameter ρ_Ph is the average density of Phobos, and a_i and a_f are the starting and final semimajor axes during any span of orbital decay. The magnitude of the surface tidal stress is proportional to the second-order tidal Love number h₂ of Phobos. Parameter h₂ can be deduced from the internal rigidity μ_i (Hurford et al. 2016), which is set to 10⁵–10⁸ Pa given the weakness of Phobos's porous interior (Andert et al. 2010). Parameter μ_o is the rigidity of the outer cohesive layer, which is set to 10⁸ Pa to be consistent with Hurford et al. (2016). The layered tidal stress model, in conjunction with the possible orbital history of Phobos, is used to reconstruct its de-orbiting stress field as Phobos spiraled closer to Mars. For simplicity, Phobos is here assumed to be a homogeneous spherical interior covered by a discrete two-layer outer regolith. Phobos's highly irregular shape and regional topography like Stickney and other craters could make the tidal stresses deviate from the symmetrical state, which has the possibility to lead to the north–south asymmetric distribution of grooves, or the mismatch could be resolved by assuming some regional compositional differences in cohesive strength. Future modeling work should include the 3D shape of Phobos to compute the surface stress field more accurately from the tidal deformation.