Discrete Element Modelling of Pit Crater Formation on Mars

Abstract

Pit craters are now recognised as being an important part of the surface morphology and structure of many planetary bodies, and are particularly remarkable on Mars. They are thought to arise from the drainage or collapse of a relatively weak surficial material into an open (or widening) void in a much stronger material below. These craters have a very distinctive expression, often presenting funnel-, cone-, or bowl-shaped geometries. Analogue models of pit crater formation produce pits that typically have steep, nearly conical cross sections, but only show the surface expression of their initiation and evolution. Numerical modelling studies of pit crater formation are limited and have produced some interesting, but nonetheless puzzling, results. Presented here is a high-resolution, 2D discrete element model of weak cover (regolith) collapse into either a static or a widening underlying void. Frictional and frictional-cohesive discrete elements are used to represent a range of probable cover rheologies. Under Martian gravitational conditions, frictional-cohesive and frictional materials both produce cone- and bowl-shaped pit craters. For a given cover thickness, the specific crater shape depends on the amount of underlying void space created for drainage. When the void space is small relative to the cover thickness, craters have bowl-shaped geometries. In contrast, when the void space is large relative to the cover thickness, craters have cone-shaped geometries with essentially planar (nearing the angle of repose) slope profiles. Frictional-cohesive materials exhibit more distinct rims than simple frictional materials and, thus, may reveal some stratigraphic layering on the pit crater walls. In an extreme case, when drainage from the overlying cover is insufficient to fill an underlying void, skylights into the deeper structure are created. This study demonstrated that pit crater walls can exhibit both angle of repose slopes and stable, gentler, collapse slopes. In addition, the simulations highlighted that pit crater depth only provides a very approximate estimate of regolith thickness. Cone-shaped pit craters gave a reasonable estimate (proxy) of regolith thickness, whereas bowl-shaped pit craters provided only a minimum estimate. Finally, it appears that fresh craters with distinct, sharp rims like those seen on Mars are only formed when the regolith had some cohesive strength. Such a weakly cohesive regolith also produced open fissures, cliffs, and faults, and exposed regolith “stratigraphy” in the uppermost part of the crater walls.

1. Introduction

Pit craters and pit crater chains are now recognised as being an important part of the surficial morphology and structure of many planetary bodies (e.g., [1,2,3,4,5,6,7]; Figure 1a,b). They have a very distinctive surface expression, often presenting funnel-, cone-, or bowl-shaped geometries, and do not possess the elevated rim, ejecta deposits, or other features that are typically associated with fresh, “young” impact craters unaffected by subsequent erosion (Figure 1a,b: see [3,4,8]). They are depressions that typically have circular-elliptical plan view shapes, and, in some cases, flat floors. They do not appear to modify the planetary surface beyond their edges. Pit crater chains are linear groups of such craters that are often, though not exclusively, associated with extensional fault traces and can occur nested within grabens (e.g., [4,8]). Exceptionally, individual craters may provide “skylights” into an underlying void ([5], Figure 1b). Their investigation is limited in natural examples as there are very few instances where we can see their subsurface structure. Thus, we typically study their surface expression and need to use analogue or numerical modelling to understand their subsurface geometry and processes. However, there are some examples where we can document natural pit craters in a more 3D sense [9,10,11,12].

Pit craters are thought to arise from the drainage or collapse of weak surficial material into an open (or evolving) void in stronger material below. The origin of such craters is somewhat contentious (see [4,7,8]). However, their morphology and geometric relations have led many authors to conclude that they are the surface traces of open voids or fractures in an underlying body due to a variety of possible formation mechanisms: extensional fracturing and dilational faulting, karstic dissolution, thermokarst processes related to sublimation or melt of ground ice, shallow dike intrusion, and lava tube collapse (e.g., [5,8,13,14,15,16,17,18]).

Analogue models of pit crater formation provide a map-view picture of their initiation and evolution, but give little insight into the internal structure of such craters (e.g., [1,2,19]). Horstman and Melosh (1989) [1] described analogue pit craters as being “a cone-shaped region with apex at the bottom the regolith and with a slope angle of about 30°”, equivalent to measured angles of repose for their regolith material, while Sims et al. (2004) [2] described them as visibly and clearly cone-shaped with steep walls. Numerical modelling studies of their formation are limited and have produced some interesting, but puzzling, results whereby the simulated pits in both weak, frictional and “strong” materials had generally convex (steepening downward) slope profiles with no distinct rim, and where the angle of repose slopes were not produced; quite unlike those observed on Earth, Mars, or in analogue modelling studies ([4]; Figure 1c). In contrast, conceptual models of pit crater drainage do not show these features, instead they typically depict cone-shaped craters with a distinct rim and angle of repose wall slopes that are planar and not downward steepening (Figure 1d). In addition, some Mars Orbiter images seem to show a kind of stratification in pit crater walls and this may indicate that, rather than being simply unconsolidated, granular materials, the cover (regolith) may well possess some cohesive strength. This hypothesis is reinforced by the occurrence of distinct, sharp rims in many pit craters and fault scarps, graben, and other features in regolith at the planetary surface.

This study presents a high-resolution (1–2 m) discrete element model of cover collapse above either a simple static or a widening underlying void (Figure 1e,f). Frictional and frictional-cohesive discrete elements were used to represent a range of probable cover rheologies. Under Martian gravitational conditions, these materials produced cone- and bowl-shaped pit craters. For a given regolith thickness, the specific pit crater geometry depends on the amount of underlying space created. No significant difference was seen between craters produced above static versus widening voids. When the void space was small relative to the cover thickness, the craters had bowl-shaped geometries. In contrast, when the void space was large relative to the cover thickness, the craters had cone-shaped geometries and crater walls had essentially planar (approaching the angle of repose) slope profiles. Frictional-cohesive materials developed craters with more distinct rims and may reveal some stratigraphic layering on pit crater walls. In extreme cases, when regolith thickness was insufficient to fill the underlying void (or was partially eroded or soluble in the case of evaporites), skylights into the deeper structure were created. Implications of these results for the interpretation of pit craters on Earth, Mars, moons, and other planets are discussed.

2. Methodology

A 2D discrete element numerical model was used to simulate the growth of (collapse) pit craters in a weak cover regolith above an underlying void with a strong basement. The modelling of the deformation to high strain was an ideal candidate for the discrete element technique, as it is well-suited to studying problems in which discontinuities (shear-zones, faults, fractures, and voids, etc.) are important. It also allowed for deformation involving unlimited relative motions of individual elements and complex, abrupt and changing boundary conditions [20,21,22,23,24,25,26], which is particularly germane to the simulation and formation of pit craters. The discrete element code used here to undertake these simulations is “cdem2D”, which was used previously to undertake a wide variety of structural geology modelling in both 2D and 3D (as cdem3D). The codes themselves were developed “in-house” and are not yet open-source, although they have been used in various external collaborations. The basic numerical methodology, and its application to a variety of problems (fault-propagation folding, orogenic wedge growth, caldera collapse, dike intrusion on Mars, viscous flow, and Gilbert deltas, etc.), were published previously [23,25,27,28,29,30,31,32]. Of particular interest here, the details of this code, applied to frictional and frictional-cohesive materials, can be found elsewhere (e.g., [24,29,32,33]). Here, the use of this modelling scheme to simulate the growth of pit craters on Mars will be illustrated.

3. Model Parameters, Set-Up and Boundary Conditions

Below, the results of two sets of simple, end-member numerical models (experiments) of pit crater formation are presented. In each case, the collapse and drainage of a (relatively) weak regolith material into a void at the regolith and basement interface (Figure 1e,f) were examined. The void was located centrally to avoid any lateral wall boundary effects. The regolith had a uniform thickness of 100 m and all models were 400 m wide. As regolith thickness on Mars is not well known, this value was chosen as a conservative minimum, and a later, thicker model will be used to examine the effect of regolith thickness. The void into which collapse occurs grows either as a result of the incremental widening of a dilational, extension fracture (Series 1 Experiments; Figure 1e), or is created instantaneously, simulating roof collapse into e.g., a lava tube (Series 2 Experiments; Figure 1f). In both experimental series, two different void geometries were considered in order to examine pit craters created, where the amount of void space relative to regolith thickness was small and large. The voids created had 2 distinct areas: with heights and widths of 100 × 40 m and 300 × 40 m, allowing for the regolith to experience limited and significant drainage. The voids that were created do not represent the actual geometry of any individual void space but, rather, create space for potential drainage of the regolith. The model represents the regolith section (the cover above the basement) and the walls and faults using approximately 19,000 discrete elements. Element radii ranged from 0.5 to 1.25 m (average radius 0.77 m) and their density was 2000 kg/m³. As the regolith on Mars is thought to include breccia and rock fragments from the local bedrock, and appears to be cohesive and stratified in places, this density was chosen as it lay between that expected for fine soil and that for intact basaltic rocks. The regolith was coloured using 12 red and yellow layers to highlight deformation, but these layers have no mechanical significance (Figure 1e,f).

The discrete elements making up the cover obey Newton’s equations of motion and interact with each other under the influence of Martian gravity. Elements experience Mohr-Coulomb frictional contact interactions with their neighbours (see [21,29] for a full description of the modelling approach). In discrete element models, such as the one used in this study, parameters such as the strength, and coefficient of friction, etc., of an assemblage are emergent properties and do not relate directly to micro-properties. They are typically assessed through the use of angle of repose and unconfined and confined biaxial numerical tests [21,25,27,34,35]. Using such tests, the frictional-cohesive discrete elements used here were found to have a bulk coefficient of friction (µ) of c. 0.70 (internal angle of friction (ø) of c. 35°) and (where present) to have a cohesion of c. 0.1 MPa. These friction and cohesion values lay within the ranges reported for natural sedimentary, granular materials and are much smaller than those typically derived from centimetre-scale laboratory samples of strong, lithified rocks (e.g., [35,36,37]), but they do lay within the ranges reported for natural unconsolidated or weakly consolidated rock masses likely similar to the Martian regolith [35,36,37].

Here, the effect of the creation of a 40 m wide vertical void in rigid basement on regolith drainage and pit crater formation is examined; this void will either grow incrementally until it reaches its final width or will appear instantaneously as a result of collapse over the roof width (Figure 1e,f). In the case of incremental void growth, this was modelled by using two rigid basement walls placed edge to edge. Dilational growth of the fracture was simulated by slowly drawing the walls 40 m apart over the first 10% of the model run, and allowing concurrent drainage of regolith into the evolving void; after the final width was achieved, simple drainage into a static void occurred. As we were dealing with frictional-cohesive materials, there was no implicit time-scale in these models. In the case of instantaneous roof collapse, the walls were placed 40 m apart at the start of the model run, the void “roof” was removed and drainage occurred immediately. The models were run for a total model simulated time of 36,000 s, sufficient to ensure that all elements reached a static equilibrium and the final configuration of the simulation was achieved.

As the experiments proceeded, the effect of the drainage of regolith elements into the underlying void was transmitted upwards into the cover, changing local element interactions and thus contact forces. As a result of this change in local forces, the internal elements were advanced to their new positions within the model by integrating their equations of motion using Newtonian physics and a velocity Verlet based numerical scheme ([38] Mora and Place, 1993). Element positions were saved during experiments to allow a detailed, high-resolution analysis of the geometry, displacement, and strain [39]. The numerical code used (cdem2D) was parallelised (using multiple computational threads) using OpenMP and was thoroughly tested and verified against the serial (only one computational thread) version [33,40]. Typical experiments like those discussed here take ~12 h to run on a DELL desktop machine with two 6-core Intel Xeon (X5650—2.66 GHz) processors allowing 24 computational threads. Model results were saved at regular intervals throughout as “time snapshots”.

For each experiment discussed herein, many different models were run with parameters similar to those described above. However, the specific experiments discussed are representative of the evolution typically observed under these boundary conditions in that they contain the reproducible, characteristic features seen in many models. The important scientific message to be taken away is not the precise location of an individual geometric feature, but rather the distinctive, reproducible behaviours that emerge from multiple experiments.

3.1. Series 1 Experiments—Progressively Opening Void (Extension Fracture)

In the first series of experiments (Figure 2), the effect of a progressively opening void (extension fracture) on an overlying regolith with either a simple frictional or a frictional-cohesive rheology is shown. As discussed above, two distinct types of void were considered—one of which had a small (potential) drainage area relative to the regolith thickness, and another which had a larger (potential) drainage area. In each case, the resultant geometries and their characteristic features at the end of the model run, when all elements were static, are examined and described.

Small relative drainage area.Figure 2a shows the final geometry as a result of the collapse of the simple frictional regolith into the smaller void and, in Figure 2b, that of the frictional-cohesive regolith. In both cases, we see that the small relative drainage area led to bowl-shaped crater geometries, sometimes exhibiting an almost flat base. None of the cases produced angle of repose (c. 35°) slopes, but the slopes of the crater walls also did not steepen downwards, in fact they flattened to a sub-horizontal central area from approximately 20–25° at the crater edges. The crater rim was less well defined in the purely frictional (granular) regolith (Figure 2a), and was very marked in the frictional-cohesive regolith (Figure 2b). Clearly, if the crater rim was ill-defined, then some appearance of steepening downward profiles in this region could be produced, but no systematic downward steepening of crater walls was observed. In addition, with frictional-cohesive regolith, we saw a rougher upper crater surface due to the formation of a number of small fault scarps. The maximum heights of these scarps were around 10 m, consistent with a low, but finite, cohesive strength. The maximum depths of the craters were of the order of approximately 30 m (compared to an initial regolith thickness of 100 m) and they had widths of approximately 150–170 m.

Large relative drainage area. Figure 2c shows the final geometry as a result of the collapse of the simple frictional regolith into the larger void, and, in Figure 2d, that of the frictional-cohesive regolith. In both cases, a larger relative drainage area led to very distinctive, cone-shaped geometries with almost planar crater walls dipping at around 30–35°. Such geometries are very different to those produced above the smaller void (compare Figure 2a,b and Figure 2c,d). The simple, frictional regolith produced quite smooth crater walls, while, in the case of the frictional-cohesive regolith, we saw a rougher upper crater surface due to small fault scarps, cliffs, and open fissures. The maximum scarp or cliff heights were around 10 m, consistent with the finite, but low, cohesive strength of the regolith material. In both cases, the crater walls met at a singular point centrally located above the infilled void. In addition, in the upper parts of the crater walls, part of the “stratigraphy” of the regolith was exposed, whereas other parts of the walls were “mantled” by drained-collapsed material lacking internal structure. The maximum depths of the craters were similar in both cases, and were of the order of around 85–90 m (compared to an initial regolith thickness of 100 m) and they had widths of approximately 230–250 m.

3.2. Series 2 Experiments—Instantaneous Roof Collapse

In the second series of experiments (Figure 3), the effect of a void, which instantaneously appeared as a result of “roof” collapse, allowing for drainage of an overlying regolith with either a simple frictional or a frictional-cohesive rheology, is shown. As before, two distinct voids were considered, one of which had a small (potential) drainage area relative to the regolith thickness, and another which had a larger (potential) area for drainage. In each case, the resultant geometries, and their characteristic features at the end of the model run when all elements are static, are examined and described.

Small relative drainage area. Figure 3a shows the result of the collapse of the simple frictional regolith into the smaller void and, in Figure 3b, that of the frictional-cohesive regolith. In both cases, as with the progressively widening void, we see that a small relative void area led to bowl-shaped geometries, with an almost flat base. Angle of repose crater wall slopes were not produced in any of the cases, and, in addition, the slopes of the crater walls did not steepen downward. The crater rim was less well-defined in the purely frictional (granular) regolith (Figure 3a), while it was very marked in the frictional-cohesive regolith (Figure 3b). In the case of the frictional-cohesive cover, we see a more stepped, rougher upper surface due to a number of small fault scarps, cliffs, and fissures. The maximum scarp heights were around 10 m, which was consistent with the low cohesive strength of the regolith. The maximum depths of the pit craters were of the order of approximately 25–30 m (compared to an initial regolith thickness of 100 m) and they had widths of approximately 135–150 m.

Large relative drainage area. Figure 3c shows the result of the collapse of the simple frictional regolith into the larger void and, in Figure 3d, that of the frictional-cohesive regolith. In both cases, as with the progressively widening void, we see that a large relative void area led to cone-shaped geometries, with almost planar sidewalls dipping at around 30–35 degrees (approaching the angle of repose). In the case of the frictional-cohesive regolith, we see that the rougher upper surface, due to small faults, cliffs, open fissures, and overhangs, was consistent with the low cohesive strength of the regolith. Such craters met at a singular point centrally located above the infilled void. This allowed us to see some stratigraphy in the upper parts of the pit craters. The maximum depths of the craters were of the order of around 85–90 m (compared to an initial regolith thickness of 100 m) and they had widths of approximately 220–230 m.

4. Discussion and Conclusions

The drainage of a (relatively) weak cover (regolith) material into an underlying void in the stronger basement, and the consequent creation of pit craters, was simulated using the discrete element method. Many of the distinctive features of pit craters seen on Mars were reproduced, including bowl-shaped geometries, cone-shaped geometries, district crater rims and angle of repose wall slopes. Whether the underlying void opens progressively (as an extension fracture) or appears instantaneously (as a result of roof collapse) made little difference to the final surficial expression of the pit crater. What made a significant difference, however, was the relative area of the void created compared to the thickness of the regolith. When the void space was small relative to cover thickness, craters had bowl-shaped geometries. In contrast, when void space was large relative to cover thickness, craters had more cone-shaped geometries with essentially planar (approaching the angle of repose) slope profiles. In every case, frictional-cohesive regolith materials exhibited more distinct rims than simple frictional materials and, thus, may reveal some stratigraphic layering (if present) in the upper pit crater walls.

Several key observations arise from these simulations. Firstly, this study emphasised that pit crater walls can exhibit both angle of repose slopes (c. 35° in this case e.g., Figure 2d) and stable, more gentle collapse slopes (c. 20° or less, Figure 3a,b). Angle of repose wall slopes tended to occur when there was significant drainage into the underlying void. Secondly, the simulations also highlighted that pit crater depth only provided a minimum estimate of regolith thickness, as there was no a priori reason to assume that void area (in 2D) would always be sufficient to drain to the regolith-basement contact. Previous studies applied geometric relationships, based on cone shapes and angle of repose slopes, to estimate pit depth and volume (e.g., [3,41]). In general, this study shows that only cone-shaped pits gave a reasonable estimate (proxy) of regolith thickness, whereas bowl-shaped pits provided, at best, a minimum estimate. Thirdly, it appears that craters with distinct, sharp rims, like those seen on Mars, were only formed when the regolith had some cohesive strength. A weakly cohesive regolith also produced open fissures, cliffs, and faults, and exposed regolith “stratigraphy” in the uppermost part of the crater walls (Figure 2b,d, Figure 3b,d and Figure 4a), while lower parts of the crater walls were mantled by disrupted, collapsed materials (Figure 4b). This is consistent with the occurrence of fault scarps and graben systems in Martian regolith ([8] Ferrill et al., 2011), and suggests that, in order to simulate and understand Martian pit craters, a frictional-cohesive cover or regolith is necessary. Finally, this study showed that bowl-, cone-, and funnel-shaped crater geometries can all exist, and that they represent a spectrum of potential geometries.

It is entirely possible that drainage from the overlying regolith may be insufficient to fill any underlying void (either the void is too large, the regolith is too thin, or the collapse material was partially or totally eroded or dissolved in the case of evaporites). In such a situation, it was suggested that skylights into the underlying void can be created [5]. These were proposed as candidate cave entrances into Martian near-surface lava tubes and volcano-tectonic fracture systems. Such Martian caves were suggested as “promising potential sites for future human habitation and astrobiology investigations” [5]. Using the same experimental set-up as the Series 2 experiments, this study showed what happens if this were to be the case for a frictional-cohesive regolith (Figure 4c). The experiment presented here was a roof collapse model, identical to the other examples, but with a very deep void (500 m) into which to drain the regolith, thus there was more volume than available regolith material. This simulation was run three times, as long as the other experiments, in order to reach a static equilibrium. From the final result we see that, indeed, a skylight into the underlying void was created under these conditions; with steep crater walls terminating at the regolith-basement interface and a partially filled void below (Figure 4c).

However, the question does arise as to why the results of these numerical simulations are so different to those presented by Smart et al. (2011) [4]. They investigated a very similar scenario to that presented here: an extension fracture underlying either a weak, unconsolidated regolith or a regolith containing strong “rock” layers (Figure 1c). They also examined a dilational fault underlying a regolith, which was not considered herein. For purely frictional, granular, weak materials, they produced some interesting, but nonetheless puzzling, results whereby the simulated pit craters had generally convex (steepening downward) slope profiles and no distinct rim; quite unlike many of those observed on Earth or on Mars (Figure 1c). Without knowing the bulk coefficient of friction of their unconsolidated material (no calibrated angle of friction for the assembly was quoted), or the bulk density of the modelled regolith, it was impossible to run comparable experiments and, thus, difficult to speculate on the origin of these differences. In order to address the idea that the Mars regolith might have some cohesive strength, they also simulated a cover which possessed a kind of mechanical stratigraphy. They did this by using several combinations of bonded elastic materials to represent moderate to very strong rock layers (cohesion ~7 MPa, and unconfined compressive strength ~22 MPa), quite unlike those suggested for Martian regolith materials. These simulations still failed to reproduce the characteristic geometries of simple Martian pit craters, apart from producing highly fractured, block-like upper crater surfaces (Figure 1c). In no case was any stratigraphic layering exposed in the crater walls. One notable difference was that the regolith used by Smart et al. was very thick (1000 m, compared to the 100 m thick regolith used here), but this scale difference should not have affected the final crater geometries in frictional-cohesive materials. To illustrate this point, a final example of a roof collapse model beneath a 1000 m thick regolith, which was run using the roof collapse experimental set-up presented herein, is presented in Figure 4d. One can clearly see that the typical crater geometries seen in the 100 m thick regolith were produced as before (compare Figure 3c and Figure 4d).

This paper is only an initial step in understanding the formation of pit craters on Mars, other planetary bodies, and their moons. It is by no means complete, however the regolith rheologies examined thus far were simple, but instructive, and reproduced many key features observed on Mars. Dilational faulting, where the void is produced as a result of a fault with a distinct throw, was also not examined. Finally, this study was undertaken in 2D, and the manner in which isolated pit craters grow, and their spacing, to form pit crater chains, coalescing in 3D to form cuspate grooves and troughs, is of great interest. These, and other topics, are the focus of current research.