Void structure stability in wet granular matter and its application to crab burrows and cometary pits

Cohesive granular matter can support stable void structures, which can universally be found in various scenes from everyday lives to space. To quantitatively characterize the stability and strength of a void structure in cohesive granular matter, we perform a simple tunnel-compression experiment with wet granular matter. In the experiment, a horizontal tunnel in a wet granular layer is vertically compressed with a slow compression rate. The experimental result suggests that the tunnel deformation can be classified into the following three types: (i) shrink, (ii) shrink with collapse, and (iii) subsidence by collapse. Using the experimental result, we estimate the stable limit of various void structures in a cohesive granular layer from crab burrows on a sandy beach to the pits observed on cometary surfaces.

On the earth, various types of natural and artificial void structures can be found in soil. Many tunnels have been constructed through mountains or under the ground to develop transportation facilities. In addition, some animals burrow into the soil to dwell in and to escape from predators. A tunnel in a sand castle constructed during children's sand play is also a type of void in soil-like material. Even in space, void structures might be ubiquitous. For example, intact lava tubes have recently been found on the Moon 1 . Moreover, pit structures observed on the surface of comets might be caused by the collapse of internal voids 2 . The mechanical properties of soils supporting these void structures are key quantities to discuss their persistence. To estimate the collapse condition, we have to understand the mechanical stability, strength, and deformation mode of tunnels. In civil engineering, the stability of tunnel structures has been studied [3][4][5][6][7] . In these studies, the stability of tunnel structures even in the excavation process has mainly been studied. In general, tunnel wall reinforcement is necessary to maintain a stable tunnel in soft soil.
Some tidal and shore animals dig unreinforced burrows without any knowledge of civil engineering. Among these, ghost crabs (Ocypode) are known to dig deep I-, J-, or Y-shaped burrows (up to ~1 m in depth) in sandy substrates 8 . In fact, previous observations revealed that unreinforced burrows can exist in unconsolidated sandy substrates 8,9 . Particularly, Seike and Nara reported that the moisture sufficiently strengthens the burrows to retain their structures 8 . However, any quantitative characterization of burrow structure strength has not been investigated.
A crab burrow in sandy shore can be regarded as a type of void structure in wet granular matter. Thus, its stability relates to the mechanical properties of wet granular matter, which usually exhibits complex behaviours [10][11][12][13][14] . Such complex behaviours of wet granular matter basically stem from the variation of the liquid structural morphology within the wet granular matter 13,14 . Although some studies exist on the mechanical characterization of wet granular matter [14][15][16][17][18][19][20][21][22][23] , little is known about the stability and/or strength of void structures formed in wet granular layers. For instance, while strength of a vertical beam of wet granular matter has been estimated 22 , the stability of a void in wet granular matter has not been investigated thus far. Therefore, in this study, we perform a simple void-compression experiment with wet granular matter. Using the experimental result, we discuss burrow stability on a sandy beach. In civil engineering, centrifugal acceleration has been used to mimic the behaviour of large-scale tunnels by small-scale laboratory experiments 3,4 . For centimetre-scale crab burrows, however, direct laboratory experiments are sufficient to directly study their mechanical properties.
The void structure in cohesive granular matter could also relate to surface terrains on comets. Since comets consist of a mixture of solid grains and ice, their surface dynamics are governed by cohesive granular mechanics. Although liquid water is basically absent in space, wet granular matter is the simplest instance of cohesive  2,24 . Although the origin of these pits is still controversial, a recent study of comet 67 P/Churyumov-Gerasimenko 2 suggests that pit structures result from the collapse of voids due to the sublimation of volatile materials in the comet. The physical mechanism governing the cometary pit size may be better understood based on our experimental results.
In the experiment, a wet granular layer with a horizontal circular tunnel is vertically compressed by a piston with a constant speed (0.5 mm s −1 ), which is sufficiently slow to neglect the inertial effect (Fig. 1). During compression, the tunnel deformation and loading force are recorded. Using these data, we characterize the mechanical property of a tunnel structure in wet granular matter. The control parameters in this experiment are the initial tunnel diameter 5 ≤ D 0 ≤ 80 (mm), grain diameter 0.1 ≤ d ≤ 0.8 (mm), initial volumetric liquid content 0.013 ≤ W 0 ≤ 0.33, and initial solid fraction (packing fraction)   φ . .

Results and Analyses
Deformation mode of a loaded void. Typical images of compressed tunnel structures are shown in Fig. 2(A-C). By the systematic experiments, we find that the deformation mode can be classified into the following three types: (i) shrink ( Fig. 2(A)), (ii) shrink with collapse ( Fig. 2(B)), and (iii) subsidence by collapse ( Fig. 2(C)). In the case of (i), a tunnel gradually shrinks without collapse. In type (ii), part of the tunnel ceiling suddenly collapses during its shrink. In type (iii), the ceiling of the tunnel completely falls down, causing a catastrophic subsidence.
From the acquired tunnel images, the area of the cross-section of the tunnel A and the corresponding diameter D (=2(A/π) 1/2 ) are calculated. In Fig. 2(D-F), temporal evolution of A is displayed. The loading force F required to compress the layer is presented in Fig. 2(G-I). In Fig. 2, the colour of the vertical dotted lines in each plot corresponds to that of the outline of the tunnel-deformation images. Since the granular layer is compressed at a constant compression rate, A monotonically decreases over time. Sharp decreases in A(t) indicate void collapse. When a tunnel shrinks without collapse, F monotonously increases. In contrast, F shows a certain maximum in the collapse cases similar to plastic deformation. Perhaps, these deformation modes might not be intuitively very surprising. Nevertheless, to the best of our knowledge, such a deformation mode classification of the compressed void in wet granular matter has not been reported in neither civil engineering nor granular physics fields.
In order to investigate which parameter significantly affects the tunnel-deformation mode, we make phase diagrams of the deformation mode as shown in Fig. 3. From these phase diagrams, one can confirm that the deformation mode is mainly dependent on d and D 0 but almost independent of W 0 and φ 0 . The large d or D 0 results in sudden collapse. Note that the deformation mode is determined by the initial diameter D 0 rather than the instantaneous diameter D. Additionally, the dimensionless factor D 0 /d cannot unify the phase diagrams. Both d and D 0 must be sufficiently small to stabilize the tunnel structure. Thus, we have to consider the individual mechanisms for the instability due to large d or D 0 . We consider these effects as follows. The d-dependent weakening of wet granular matter is not surprising. Because the larger d yields the larger radius of curvature of capillary bridges among grains, it weakens the cohesive effect due to the decrease of Laplace pressure. However, D 0 -dependent weakening of tunnel structure cannot be easily understood. The simplest way to characterize the D 0 dependence could be a balance between a mass-defect-based stress and material strength. Here, the mass-defect-based stress Δσ 0 can be estimated as Δσ 0 = ρgV void /S proj , where ρ, g, V void , and S porj correspond to the bulk density of the wet granular layer, gravitational acceleration, void volume, and area of the void projected to a horizontal (perpendicular to gravity) surface, respectively. If Δσ 0 , which roughly corresponds to the defect of static soil pressure ( σ ρ Δ~gD 0 0 ), exceeds the material strength, the void structure becomes unstable causing the collapse even with very small additional loading. Namely, Δσ 0 effectively corresponds to the initial loading stress due to void existence. If there is not any void structure, the initial loading is absent; Δσ 0 = 0. When Δσ 0 exceeds the material strength, the structure becomes unstable, leading to collapse. To evaluate this condition, we have to estimate the effective material strength of a void in wet granular matter.

Strength of tunnel structures.
To estimate the tunnel-based strength, we employ a simple model for the maximum shear stress applied to a tunnel structure in soil 25 . In this model, the maximum shear stress acting on the tunnel structure, τ, is obtained from a static force balance in the vertical direction as 25 where σ ex = F/XY is the stress computed from the external loading force F acting on the top surface area XY (Fig. 1). C is thickness of cover upon a tunnel (Fig. 4(A)). As shown in Fig. 4(A), the two-dimensional cross-section of the tunnel is considered in the model, and τ corresponds to the shear stress acting on the plane of 45° tilted from the vertical direction. In this model, we simply consider the magnitude of stress τ at the top part of tunnel. We assume that when τ exceeds the material strength, deformation or fracturing is induced. Figure 4(B) shows τ computed using Eq. (1) as a function of D/D 0 . The tunnel diameter D is measured from the acquired images. The bulk density ρ is computed from the volume of the compressed wet granular layer. Other parameters are set by initial condition or constant. At the beginning of loading (D/D 0 ≃ 1), τ increases rapidly. When D 0 is small (21 mm), continuous shrinking of the tunnel occurs and τ reaches ≃ 3 kPa. In the subsidence case (D 0 = 75 mm), however, τ shows a sharp peak due to the collapse without significant shrink. The maximum value of τ in this case is less than 2 kPa. The shrink-and-collapse phase (D 0 = 40 mm) shows intermediate behaviour between them. To estimate the stability limit, we define the characteristic strength by the stress value at which the tunnel deformation begins, τ yield . We define the onset of deformation by the threshold of the diminution rate, |dA/dt| ≥ 1 mm 2 s −1 . An example of dA/dt is shown in Fig. 4(C) in which the sudden contraction of A can clearly be confirmed above the threshold level. Since most of the dA/dt data show a similar trend, we choose this specific threshold value. However, the τ yield value depends on the threshold value. This arbitrariness actually results in the large variance of τ yield as discussed later. In addition, we define another strength of the tunnel structure τ max by the peak of τ. Peaks in stress curves (Fig. 4(B)) for the shrink-with-collapse and subsidence-by-collapse phases come from the actual failure of the tunnel shape, as shown in Fig. 2(H,I). However, the shrink phase does not show any sudden collapse and corresponding decrease in F (Fig. 2(G)). The peak of the stress curve for the shrink-phase tunnel shown in Fig. 4(B) comes from the decrease in D. In Eq. (1), the smaller D results in smaller τ. When this effect overcomes the effect of σ ex increase, τ decreases. In other words, to keep τ sufficiently small (below the strength level), the tunnel shrinks in the shrink phase. Therefore, the maximum of τ in the shrink phase can also represent the bulk strength of wet granular matter. In the following analysis, we mainly focus on the shrink-phase tunnels to estimate the strengths using Eq. (1). The dependencies of τ yield on W 0 and φ 0 are shown in Fig. 5(A,B), respectively. τ yield distributes around 0.5 kPa. The order of these experimentally obtained strength values (~10 0 kPa) is consistent with previous studies on strength of wet granular matter 11 . τ yield increases with W 0 at small W 0 regime and shows a peak around W 0 = 0.15. Liquid bridges among grains strengthens the granular layer in W 0 ≤ 0.15. Then, τ yield decreases at a larger W 0 regime because liquid bridges are replaced by liquid clusters, which provide a lubrication effect 20 . However, W 0 dependence of τ yield is very weak; thus, we can assume it is approximately a constant. τ yield does not clearly depend on φ 0 . The large variance of τ yield vs. φ 0 (Fig. 5 (B)) could originate from the effect of the fixed threshold for computing τ yield . Although the value of τ yield depends on the threshold of dA/dt, the qualitative behaviour and the order of magnitude of τ yield are almost independent of the threshold value. To evaluate the instability of the void structure, the comparison between Δσ 0 and these strength values would be helpful. By substituting typical values ρ = 1.4 × 10 3 kg m −3 and V void /S proj = πD 0 /4 with the typical initial diameter of collapsing tunnel D 0 = 50 mm into the form of Δσ 0 , we obtain σ Δ .  0 54 0 kPa. This value is actually close to τ yield ≃ 0.5 kPa. That is, when a dimensionless number χ = Δσ 0 /τ yield exceeds unity, the initial void structure is unstable and subsidence by collapse could be induced by small amount of additional loading. Since τ yield is relatively independent of W 0 , D 0 is the principal variable in χ. That is the reason why the phase boundaries in Fig. 3 are almost horizontal. When χ is sufficiently less than unity, on the other hand, continuous shrink is induced by the compression. Figure 5 (C,D) show the variations of τ max as a function of W 0 and φ 0 . The value of τ max is approximately in the range of 1-10 kPa. τ max monotonically increases with W 0 . This trend differs from that of τ yield because the effect of compaction is not negligible for τ max . When the granular layer is loaded, pressed grains contribute to not only filling a tunnel hole but also increasing packing fraction of the matrix (bulk) layer. Since the wetting liquid allows for grains to move smoothly, the granular layer can be easily compacted. Whereas this matrix-compression effect is limited in the early stage of deformation, it plays a certain role in the late stage determining τ max . As shown in Fig. 5(D), τ max increases with φ 0 . This trend likely originates from the increase of contact points among grains where new liquid bridges can be formed. Some empirical fitting forms are also presented in Fig. 5 only to guide the eye. Since the variations of D 0 and d are limited in this experiment, D 0 and d dependencies of τ yield and τ max can only be roughly discussed. While the data are not shown here, we find that both τ yield and τ max are almost independent of D 0 and are decreasing functions of d. Details of the systematic analyses of the entire data set and the empirical forms for estimating strength values are found in 26 .

Discussion
Using the experimental result, the upper limit of crab burrow size can be discussed. We have also investigated crab burrows at a sandy beach in Tsu, Mie prefecture, Japan. The average diameter of the actual crab burrows we found is 26.4 mm, and the standard deviation is 6.2 mm. In addition, the water content W of the region at which crab burrows exist is in the range of 0.1 ≤ W ≤ 0.4. These observations are consistent with those of previous studies 8,9 . From the observational results and phase diagrams in Fig. 3, we can conclude that burrows are safety in terms of collapse prevention. Put differently, ghost crabs make sufficiently small burrows (D 0 < 50 mm) to reduce the risk of sudden collapse. Since the phase boundary does not significantly depend on W 0 (at least in the range of W at the actual sandy beach, 0.1 ≤ W ≤ 0.4 (Fig. 5(A))), the water content is not very important as long as the layer is somewhat wet. This W 0 -insensitivity is advantageous for crabs because the water content on beach could vary depending on the weather/tide 27 . The results presented thus far could depend on the grain shape. According to the preliminary results of the experiment with sand grains, the strength of sand qualitatively shows the similar tendency to the current experimental result obtained by glass beads. The structure of the phase diagram shown in Fig. 3 is almost identical in the sand case. A much more detailed comparison between the experimental result and actual crab burrows on the basis of field survey will be presented elsewhere.
By considering the buckling instability, the strength of a wet granular pillar was estimated 22 . According to their result, the typical strength of a wet granular pillar is on the order of MPa. This value is much greater than the strength estimated in this study. This difference is rather natural since the dense pillar structure is much stronger than the void structure. The strength obtained in this study must be used to properly discuss the stability of sand tunnel structures.
The current experimental result also enables us to discuss the size of pits observed on comet 67 P. Vincent et al. reported that the size of pits ranges from 50 m to 300 m 2 . In that study, they only discussed the relation between critical ceiling thickness and void diameter by assuming a simple sinkhole model. However, in their model, the size range is arbitrary. In our experiment, it turns out that the small void is stable. Moreover, the small voids result in the continuous shrink without collapse. Namely, the small voids cannot produce pits. Thus, the lower limit of size of pits could be estimated by the comparison between Δσ 0 and the strength of materials. According to Vincent et al. 2 , the typical values for comet 67 P are ρ = 470 kg m −3 , g = 5 × 10 −4 m s −2 , and tensile strength of surface material τ com ≃ 50 Pa. Here, we consider the collapse condition by calculating stress balance among the effective loading by the void structure Δσ 0 , the yield stress τ yield , and the maximum material strength τ max with a help of the collapse condition χ = Δσ 0 /τ yield ≃ 1. By assuming a spherical void shape with the lower-limit diameter of collapsing void, D 0 = 50 m, Δσ 0 becomes Pa. To consistently explain the pit subsidence by void collapse, the effective τ yield on the surface of a comet should roughly be identical to Δσ 0 = 8 Pa from the criterion χ ≃ 1. From the experimental result, the ratio between material strength τ max and yield strength τ yield is τ max /τ yield ≃ 3-6 ( Fig. 4(B) and τ yield ≃ 0.5 kPa). This value is approximately consistent with the observation, τ com /τ yield = 50/8 ≃ 6 by assuming that τ max is governed by tensile strength. Namely, the lower limit of pit size could be determined by the stability limit of the void size in the cohesive granular layer.

Conclusion
In summary, we conducted a simple experiment of shrink and collapse of voids in compressed wet granular matter. The phase boundary between shrink and collapse mainly depends on grain size d and initial void size D 0 . Particularly, the tunnel structure with larger d or D 0 tends to collapse. From the experimental result, we revealed that the crab burrows are sufficiently small to prevent the subsidence hazard, and the lower-limit size of pits on comet 67 P could be governed by a similar collapse mechanism.

Methods
In the experiment, we observe how the tunnel structure in the wet granular layer deforms when it is uniformly compressed from the top of the layer. The wet granular matter is prepared by rotating a mixture of grains and liquid (tap water, ρ liquid = 1000 kg m −3 ) for 10 minutes with 100 rpm on the pot mill rotation table (Nittokagaku, ANZ). The tunnel structure is made by extracting a solid cylinder (of diameter D 0 ) from the wet granular layer packed in an acrylic container with holes. Note that the solid cylinder is placed before preparing the wet granular layer and withdrawn after packing. To neglect the aging effect such as draining or evaporation of water, the granular layer is quickly prepared, and the prepared granular layer is compressed right after the preparation. We confirmed that the evaporation is sufficiently slow. Regarding draining, its effect is negligible due to the capillary force when the liquid content is not large. However, in the large liquid-content case, draining might affect the result. To minimize that effect, we quickly perform the experiment.
In this experiment, we control the following parameters: initial liquid content W 0 , initial solid fraction φ 0 , initial diameter of the tunnel D 0 , and grain size d. W 0 is defined as the volume of liquid V liquid divided by the total volume of the wet granular layer including solid, liquid, and void parts, V total ; W 0 = V liquid /V total . On the other hand, φ 0 is defined as φ 0 = V grain /V total , where V grain is the volume occupied by (dry) grains. The experiments are carried out 3 times with identical initial conditions. We perform experiments over a range of W 0 from 0.013 to 0.33 and φ 0 from 0.44 to 0.56. Various size containers to vary D 0 are used: (D 0 , X, Y, Z) = (5 mm, 75 mm, 40 mm, 100 mm), (10 mm, 80 mm, 40 mm, 105 mm), (20 mm, 100 mm, 40 mm, 150 mm), (40 mm, 110 mm, 40 mm, 135 mm), and (80 mm, 150 mm, 40 mm, 175 mm). See Fig. 1 for the definition of dimensions X, Y, and Z. The ceiling thickness of the granular layer above the tunnel C is fixed at 20.6 ± 2.0 mm at the initial state. Although there might be a certain wall effect in the experiment, we fix the depth of the cell (Y = 40 mm) to minimize the effect of wall effect variation. Moreover, Y = 40 mm is sufficiently thick to neglect the frictional wall support (so-called Janssen effect, see, e.g., Sec. 3.7 in 28 ) since the layer thickness above the tunnel (C = 20 mm) is smaller than Y. Grains used in this experiment are glass beads (true density: ρ grain = 2500 kg/m 3 , grain size: d = 0.1, 0.4, or 0.8 mm, AS-ONE). During the compression, movies of the cross-section of deforming tunnel is acquired by a CMOS camera (Sentech, STC-MCCM401U3V (4 M Color)) at a rate of 2 frames per second by using transmitted light. We also measure the force F to compress the wet granular layer and displacement of the loading top plate using a universal testing machine (Shimadzu, AG-X). In this experiment, loading rate of the top plate is fixed at 0.5 mm s −1 . We checked that this loading rate is sufficiently slow to neglect the inertial effect. The behaviour does not depend on a loading rate below the speed of 0.5 mm s −1 . This loading rate is the fastest one in the quasi-static regime. To avoid the aging effect, the faster loading rate is better. Therefore, we employ this rate.