Experiments on seepage-triggered cliff landslides using cohesive wet sand – Progress in Earth and Planetary Science

Redistribution of effective stress with seepage and the resulting failure

Experiments show that for tall ((H / L ge 0.89)) cliffs, the (sigma _{mathrm{z:back}}) (near the back gate) rises whereas the (sigma _{mathrm{z:front}}) (near the cliff toe) falls with seepage (Fig. 9a), and that (sigma _{mathrm{z:front}}) becomes sub-lithostatic at (t = T_{mathrm{seep}}) (Fig. 9b). The rise in (sigma _{mathrm{z:back}}) cannot be explained by the rise in u alone. The fall in (sigma _{mathrm{z:front}}) starts prior to the rise in u, measured at the same horizontal distance from the hydraulic head, indicating that the groundwater has not yet arrived. These observations indicate a redistribution of effective stress (sigma ‘), arising from a change in arching, that is related to the loss of cohesion with seepage and the development of the slip plane.

Figure 12a is a schematic diagram showing the stress change inferred from the measurements. The figure models the situation shown in Fig. 6 with (H = 20) cm and (h_{mathrm{w}} = 8) cm. At (t = 0) s, the (sigma _{mathrm{z:front}}) is super-lithostatic whereas the (sigma _{mathrm{z:back}}) is sub-lithostatic, indicating an arching. At (t sim T_{mathrm{slide}}), we observe that tension cracks formed near the cliff top, indicating that an active failure by horizontal extension occurred. Tension cracks form because cohesive wet sand resists tensile stress (sigma) ((< 0)) estimated to be (sigma > – c / mu sim – 2300) Pa ((S = 0.3)) (Andreotti et al. 2013). Such cracks only form near the cliff head, above the phreatic surface, as we showed from Eq. (3). On the other hand, near the cliff toe, we infer that the sand is compressed horizontally by the pressure of the hydraulic head, such that a passive failure occurred. This stress pattern is consistent with that inferred from field landslides, where extension and compression features form at the head and toe, respectively (Whitely et al. 2019), and a slip plane which is steeper at the head (Carson 1971; Matsukura 2008); a consequence of steeper failure angle (measured from the horizontal) under tension, compared to that under compression (Fig. 12a, right).

a A schematic diagram to scale showing the inferred stress changes with seepage. The long axis of a rectangle indicates the major principal stress (sigma), defined positive when the arrow is directed inward. The figure models the situation of run 14 ((H = 20.0) cm, (L = 13.5) cm, (h_{mathrm{w}} = 8) cm shown in Fig. 6). A red square (circle) indicates the pressure sensor attached to the base (side wall) and measures the (sigma _{mathrm{z}}) ((sigma _{mathrm{y}})). Left: at the beginning of seepage ((t = 0) s). Right: at (t = T_{mathrm{slide}}). Near the cliff head, (sigma _{mathrm{z}}) is the inferred major principal stress, such that active failure (horizontal tension, (sigma _{mathrm{x}} < 0)) occurs and a tension crack forms. At the cliff toe, the (sigma _mathrm{x}) is the inferred major principal stress, such that passive failure (horizontal compression, (sigma _{mathrm{x}} > 0)) occurs. The failure angle is steeper at the head (see text). b Effective Mohr’s circle representation of a passive failure at the cliff toe shown in a right. A black line is the Mohr–Coulomb criterion for saturated ((S = 1)) sand (see Additional file 4: Table S1 for (mu ‘), (c’)). (b simeq 25^{circ }) is the failure angle measured from the horizontal for the experimental case

The pattern of (sigma ‘) in granular matter is commonly described using force chains or networks (Andreotti et al. 2013; De Blasio 2011). Figure 9a suggests that at (t = 0) s, near the cliff toe the force chains are directed vertical, whereas near the back gate they are directed toward the back gate. Such enhanced wall-effect for unsaturated granular matter, compared to saturated granular matter, has been shown by Nowak et al. (2005). When seepage starts from the rear of the cliff, the cohesion below the phreatic surface drops which we consider weakens the force chains directed toward the back gate. We infer that the force chain redirection toward the base contributed to the rise in (sigma _{mathrm{z:back}}). The fall in (sigma _{mathrm{z:front}}) occurs simultaneously with the rise in (sigma _{mathrm{z:back}}), from which we infer that this is also a result of the force chain redirection. A possible explanation is that the force chains near the cliff toe redirected toward the side walls (in the y-direction), as a consequence of horizontal compression (in the x-direction, see Fig. 12a right). Such redirection toward the wall, perpendicular to the direction of imposed stress, is a property of granular matter as imaged by photoelasticity (Duran 2000). We infer that this redirection supported the weight of the wet sand and thereby reduced the (sigma _{mathrm{z}}). Indeed a complementary increase of (sigma _{mathrm{y}}) measured at the side wall (at a height 3.5 cm) is evident in stages I–II of Fig. 6b. Here the (sigma _{mathrm{y}}) rises with seepage, exceeding the expected rise from the u. We also note that the (sigma _{mathrm{z:front}}) recovers and rises in stage III as evident in Figs. 5b, 6b and 7b. We infer that this stress recovery is a consequence of the cliff detaching from the side walls which weakened the force chains directed in the y-direction.

A consequence of the fall in (sigma _{mathrm{z:front}}) is that it promotes failure near the cliff toe, which we explain below. Figure 12b shows the effective Mohr’s circle representation of the passive failure condition in the (x – z) plane near the cliff toe at (t = T_{mathrm{slide}}) of Fig. 12a. Here we compare the following 2 cases. The first is when the (sigma _{mathrm{z}}) is at a lithostatic value. From (sigma _0 simeq 3109) Pa [Eq. (6)] and (u simeq 784) Pa [Eq. (7)], we obtain the (sigma _{mathrm{z}}’ simeq 2325) Pa [Eq. (4)]. From the Mohr’s circle, we obtain (sigma _{mathrm{x}}’ simeq 11616) Pa (coefficient of passive earth pressure (K_{mathrm{p}} = sigma _{mathrm{x}}’ / sigma _{mathrm{z}}’ simeq 5.0)) for a failure to occur. The second is when we use the measured (sigma _{mathrm{z}} simeq 953) Pa ((< sigma _0)), and (u simeq 688) Pa at (t = T_{mathrm{slide}} simeq 308) s, from which we obtain (sigma _{mathrm{z}}’ simeq 265) Pa. From the Mohr’s circle, we obtain (sigma _{mathrm{x}}’ simeq 1838) Pa ((K_{mathrm{p}} simeq 6.9)) for a failure to occur. Comparing the (sigma _{mathrm{x}}’) values for these 2 cases, we find that the (sigma _{mathrm{x}}’) for the experimental case is reduced to (sim 16%) of the lithostatic case. This shows that the fall in (sigma _{mathrm{z:front}}) promotes failure.

A permeable flow model for seepage

Here we model seepage and calculate the (T_{mathrm{seep}}) and compare it with the measured (T_{mathrm{seep}}) shown in Fig. 10c, d. We assume a permeable flow driven by a horizontal pressure gradient (the Dupuit approximation), whose velocity is determined by Darcy’s law (Bear 1972; Turcotte and Schubert 2014). We use this simple model which assumes 1-D (horizontal) flow because our goal is to clarify whether a simple model can explain the measured (T_{mathrm{seep}}) and its (h_{mathrm{w}})-dependence. The nonlinear diffusion equation, which governs the height h(t, x) of the phreatic surface, is expressed as

where k is the permeability, (eta = 1.0 times 10^{-3}) Pa s is the viscosity of water at (20,^{circ })C, and x is the horizontal coordinate measured from the origin ((x = z = 0)) defined at the cliff toe. For k we use the Kozeny–Carman formula (Bear 1972; Mavko et al. 2019),

where K is an empirical constant that depends on the particle shape and tortuosity (the ratio of the total flow-path length to the sample length). Previous works have suggested (K = 180) for spherical (Bear 1972), and (K = 4000) for aspherical particles (McKenzie 1984). For the (phi) and d we use (phi = 0.51) (Table 2) and (d = 0.2) mm. The initial condition is (h = 0) cm in the cliff ((0< x < 13.5) cm) and the boundary condition is (h(t, x = 13.5) cm) (= h_{mathrm{w}}) at the constant hydraulic head. We numerically integrate Eq. (10) until the water seeps out from the cliff toe. The seepage time (T_{mathrm{seep}}) is defined when (h(t, x = 0) (text{ cm})) at the cliff toe exceeds the d.

First we consider the (h_{mathrm{w}}) dependence shown in Fig. 10c. We calculated the (T_{mathrm{seep}}) by varying the (h_{mathrm{w}}) using 2 values of k, (k_1 = 5.0 times 10^{-11}) (text{ m}^2) and (k_2 = 1.8 times 10^{-11}) (text{ m}^2), so that the calculated (T_{mathrm{seep}}) (cyan curves in Fig. 10c) bound the measured (T_{mathrm{seep}}) (blue (bigtriangleup)). The curves explain the (h_{mathrm{w}}) dependence well (see also Additional file 4: Fig. S5 for a semi-log version). To evaluate whether these values for k are reasonable, we calculated k using Eq. (11) by substituting (phi = 0.51), (d = 0.2) mm, and using (K = 180) and (K = 4000) as cited above, and obtained (k = 1.0 times 10^{-10}) (text{ m}^2) and (k = 4.5 times 10^{-12}) (text{ m}^2), respectively. This range of k overlaps (k_1) and (k_2), indicating that their values are physically reasonable.

Next we consider the H dependence shown in Fig. 10d. The (T_{mathrm{seep}}) increasing with H suggests that k decreased with H. The k depends on both (phi) and K [Eq. (11)]. The initial (phi) does not depend on H (Table 2), and although the (phi) increases by 0.15–0.75% through compaction during seepage (Additional file 4: Fig. S6), no systematic H dependence is found. Here we note that effective stress (bar{sigma ‘_{mathrm{z}}}), which is averaged horizontally (using (sigma ‘_{mathrm{z:front}}), (sigma ‘_{mathrm{z:back}})), and over a time span of (0< t < T_{mathrm{seep}}), increases with H (Additional file 4: Fig. S7), as expected from lithostacy. Thus, we infer that the total flow-path length (the tortuosity, hence K) increased with (bar{sigma ‘_{mathrm{z}}}) (hence H), because the sand particles were in a closer contact, and that this resulted in the k to decrease.

Darcy’s law, used in Eq. (10), assumes a laminar flow (Bear 1972; Turcotte and Schubert 2014). Here we evaluate whether this assumption holds by calculating the Reynolds number (Re = rho _{mathrm{w}} {bar{v}}_{mathrm{x}} d / eta), where ({bar{v}}_{mathrm{x}}) is the time-averaged seepage velocity. Using an upper bound of ({bar{v}}_{mathrm{x}} sim 3) mm/s, which is the fastest case with (h_{mathrm{w}} = 15) cm (run 7), we obtain an upper bound estimate of (Re sim 0.6). A value of (Re < {mathcal {O}}{(1)}) validates assuming a laminar flow (Bear 1972) for the seepage in our experiments.

We note, however, that the flow assumed here is simplified because it neglects the vertical flow. The flow becomes 2-D when when dh/dx becomes of the order of unity (Turcotte and Schubert 2014), which is realized at the front of the phreatic surface. Vertical flow is also driven by suction. Here the resulting capillary rise (h_{mathrm{c}}) scales as (h_{mathrm{c}} propto gamma / (1 – phi ) rho _{mathrm{w}} g d), and using an empirical formula in Terzaghi et al. (1996), we estimate to be around (h_{mathrm{c}} sim 10) cm. A fair agreement with the measurements, despite neglecting the vertical flow, may be because the (T_{mathrm{seep}}) was measured at the cliff toe, just above the impermeable base, where the vertical velocity approaches zero.

Comparison with stability analysis

Here we evaluate the cliff stability and estimate the geometry of the most unstable slip plane, and compare them with the experiments. Same as in Sect. 6.2, we use a simple model to clarify the extent of explaining the experiments, which we consider is meaningful before seeking more complex models. Such analysis in conjunction with the experiments has been conducted previously (Fox et al. 2007; Tohari et al. 2007) which we compare with our result at the end of this section.

We use the standard method of slices (Bodó and Jones 2013; De Blasio 2011; Terzaghi et al. 1996) as shown schematically in Fig. 13. We assume a circular slip plane passing through the cliff toe and a parabolic phreatic surface h(x) (Dupuit parabola), which is a solution for a steady flow with a boundary condition (h(x = 0) cm) = 0 cm (Bear 1972; Turcotte and Schubert 2014). A rationale for assuming a circular slip, which is steeper at the head, is because it is consistent with the horizontal tension (compression) at the cliff head (toe) of the landslides (Sect. 6.1). A factor of safety (F_{mathrm{s}}) is defined as (F_{mathrm{s}} = M_{mathrm{R}} / M_{mathrm{D}}), where (M_{mathrm{R}}) is the resisting moment and (M_{mathrm{D}}) is the disturbing moment for a given slip plane. Thus, the cliff is unstable (stable) when (F_{mathrm{s}} < 1) ((F_{mathrm{s}} > 1)).

A schematic diagram showing the geometry of the method of slices used to evaluate cliff stability. The slip plane is circular (center C, radius R), passing through the origin (O). The landslide head length (L_{mathrm{head}}) is measured from the intercept of the slip circle at the top or rear plane of the cliff. The phreatic surface h(x) is parabolic, such that (h(0) = 0) and (h(L) = h_{mathrm{w}}). Wet sand is saturated (unsaturated) below (above) the h(x), and corresponding values of (rho), (mu ‘), (c’) are used. (Delta x) ((= 1.35) mm) is the slice width. (z_{mathrm{i}}), (l_mathrm{i}), and (delta _{mathrm{i}}) are the height, length, and slope angle of the slip plane, respectively, at the slice number i. (s_{mathrm{i}}) and (tau _{mathrm{i}}) are the resisting and disturbing tangential stress, respectively, at the slip plane

(M_{mathrm{R}}) and (M_{mathrm{D}}) are evaluated as follows. The cliff is divided into vertical slices, each with a width of (Delta x). The height, length, and slope angle of the slip plane at each slice are defined as (z_i), (l_i), and (delta _i), respectively, where i is the slice number. The resisting and disturbing tangential stresses at each slice are defined as (s_i) and (tau _i), respectively. (F_{mathrm{s}}) is then expressed as (Bodó and Jones 2013),

Here (w_i = rho _i g (H – z_i) Delta x) is the weight of wet sand (H: cliff height) evaluated at height (z_i), and (u_i) is the pore water pressure. We assume that the wet sand is unsaturated (saturated) above (below) the phreatic surface (height (h_i)), such that the (u_i) is zero (non-zero: (u_i = rho _{mathrm{w}} g (h_i – z_i))). Thus, we use (rho _i = 1488) (text{ kg/m}^3), (mu _i’ = 0.60), and (c_i’ = 1362) Pa ((S = 0.3) case) above the (h_i), and (rho _i = 1821) (text{ kg/m}^3), (mu _i’ = 0.86), and (c_i’ = 133) Pa below the (h_i) (Additional file 4: Table S1). We vary the center C of the slip circle (and the circle radius R), and determine the C which minimizes the (F_{mathrm{s}}). An example of such calculation is shown in Fig. 14a. From the slip circle, we determine the landslide head length (L_{mathrm{head}} (le L)) as shown in Fig. 13. The calculation is repeated for different pairs of (h_{mathrm{w}}) ((le H)) and H to study how the (F_{mathrm{s}}) and (L_{mathrm{head}}) depend on these 2 parameters.

a An example of a factor of safety (F_{mathrm{s}}) calculation in which the center of slip circle (C) is varied. A black rectangle indicates the cliff ((H = 20) cm, (L = 13.5) cm). A blue curve indicates the phreatic surface with a hydraulic head (h_{mathrm{w}} = 12) cm at (x = 13.5) cm. Filled circles indicate the C, and their colors represent the calculated (F_{mathrm{s}}) (see legend for their values). A red circle indicates the C at which the (F_{mathrm{s}}) is the minimum ((F_{mathrm{s}} = 0.99)). Broken red lines are the slip circle radii, and a red arc is the slip circle sector. b (F_{mathrm{s}}) versus (h_{mathrm{w}}) calculated for cliffs with 2 different H values (see legend). For each H, the (h_{mathrm{w}}) is varied in the range of (h_{mathrm{w}} le H). A red (bigcirc) corresponds to the minimum (F_{mathrm{s}}) in a

Now we compare the calculations with the experiments. First we consider the (F_{mathrm{s}}). Figure 14b shows that (F_{mathrm{s}}) decreases with (h_{mathrm{w}}), which is consistent with the existence of a critical (h_{mathrm{w}}) above which a slide is triggered. Here we demonstrate that the order of magnitude drop in c upon saturation is important. For the (h_{mathrm{w}} = 12) cm, (H = 20) cm case shown in Fig. 14a, if we assume that the c remains unchanged at its unsaturated value below the phreatic surface, then we obtain (F_{mathrm{s}} = 2.14), a factor of (sim 2) larger than the (F_{mathrm{s}} = 0.99) (red (bigcirc) in Fig. 14b) obtained, which included the drop in c. The calculated critical (h_{mathrm{w}} simeq 12) cm for a slide to be triggered, however, is overestimated compared to the critical (h_{mathrm{w}} sim 2) cm of the experiments. We consider that there are 2 main reasons for this discrepancy, both of which overestimate the (M_{mathrm{R}}) in Eq. (12). First is the neglect of capillary rise, estimated to be around (h_{mathrm{c}} sim 10) cm (Sect. 6.2), which causes the cohesion (c_i) of the wet sand above the assumed phreatic surface, to fall. Second is the neglect of (sigma _{mathrm{z:front}}) becoming sub-lithostatic near the cliff toe at (t sim T_{mathrm{seep}}), which promotes failure (Fig. 12b). In the analysis, lithostacy was assumed, which overestimates the (w_i) (hence frictional resistance) near the cliff toe where the (delta _i) is small.

Next we consider the (L_{mathrm{head}}). In Fig. 8b, we indicated the (L_{mathrm{head}}) calculated by varying the H for constant (h_{mathrm{w}} = 2) cm and (L = 13.5) cm. The figure shows that the calculated (L_{mathrm{head}}) increases with H (i.e., the slide becomes deep seated), in accord with the experiments. We also calculated the (L_{mathrm{head}}) by varying the (h_{mathrm{w}}) at a constant (H = 20) cm, and find (L_{mathrm{head}} = 13.5) cm for (h_{mathrm{w}} le 18) cm, indicating that the (L_{mathrm{head}}) depends little on the (h_{mathrm{w}}), which also agrees with the experiments. Furthermore, we compared the measured and calculated (L_{mathrm{head}}) for the (L = 8.5) cm case (run 18) and found (L_{mathrm{head}} = 8.5) cm for both. Our analysis thus shows that (F_{mathrm{s}}) is overestimated, whereas the (L_{mathrm{head}}) increasing with H agrees with the experiments.

Here we also interpret the result shown in Fig. 10c (also Additional file 4: Fig. S5), where the time delay (T_{mathrm{slide}} – T_{mathrm{seep}}) (i.e., the time span of stage II: (T_{mathrm{seep}}< t < T_{mathrm{slide}})) becomes longer when the (h_{mathrm{w}}) is lower. During stage II, we observe that the sand compacts (Figs. 5b, 6b). We consider that this is a consequence of the sand becoming increasingly saturated. When the wet sand becomes saturated ((S = 1)), it loses capillary cohesion (Fig. 1d), and compacts under its own weight, which was verified from rheometry conducted under an imposed normal stress (see legend of Additional file 4: Table S1). When the (h_{mathrm{w}}) is low, the cliff becomes stable (Fig. 3b, consistent with large (F_{mathrm{s}}) at small (h_{mathrm{w}}) in Fig. 14b). A longer time delay can be interpreted as a result of extra time needed for the wet sand to become saturated, not only below the phreatic surface but also above it by capillary rise (so that the (F_{mathrm{s}}) decreases), to trigger a slide. On the other hand when the (h_{mathrm{w}}) is high, the cliff becomes unstable. A slide can be triggered at (t sim T_{mathrm{seep}}), without a further increase in S.

Finally, we compare with previous works. Tohari et al. (2007) conducted landslide experiments triggered by rainfall and groundwater rise and showed that the (F_{mathrm{s}}) is overestimated unless the loss of strength by saturation below the tension crack was included, a result similar to ours. Fox et al. (2007) conducted experiments in which groundwater ((h_{mathrm{w}} / H = 0.86)) seeps into a quasi 2-D bank (cliff and slopes), consisting of sand. They found that the seepage erodes the bank toe and forms an undercutting for steep slopes, after which they collapse, a result broadly similar to our experiments, but different in detail. They compared the slip geometry of the experiments with that obtained from the stability analysis and showed that the analysis could not predict the failure geometry at small slope angles because the slip becomes non-circular.

Implications for field landslides and their caveats

Now we discuss the implications for field landslides from the perspective of the 3 main issues which we addressed in this work. First, we consider the loss of cohesion c, by an order of magnitude, with seepage. The c and (H_{mathrm{max}}) [Eq. (3)] of soil consisting of silt and clay, which are finer compared to the sand we use, also drops by an order of magnitude with S (Sect. 4). Thus, we consider that our experiments are applicable to field landslides whose soil consists of fine particles, and those which are shallow and steep, or at low-gravity (Sect. 1). We note that wet snow is similarly cohesive (Bartelt et al. 2015) and loses cohesion by an order of magnitude when its water content increases (Kamiisi et al. 2009). Thus, snow avalanches triggered by seepage or melting may also share similarities with our experiments.

Second we consider the measurements of total stress (sigma), which are less common than the measurements of pore water pressure u. We found that the initial (sigma _{mathrm{z}}) of a tall cliff deviates from lithostacy (Fig. 9a), becoming high near the cliff toe and low near the back gate, and that lithostatic assumption is one reason for the factor of safety (F_{mathrm{s}}) being overestimated (Fig. 14b). (sigma _{mathrm{z}}) has been measured in the field (Watts and Charles 1988) but less commonly compared to the horizontal (sigma) (Araiba and Suemine 1998; Askarinejad et al. 2018; Brackely and Sanders 1992; Cislaghi et al. 2019). Our experiments show that for cohesive soils, one should be cautious about assuming lithostacy and therefore it is meaningful to measure (sigma _{mathrm{z}}). We showed that (sigma _{mathrm{z}}) monitoring detects stress redistribution before the groundwater arrives. This suggests that the (sigma _{mathrm{z}}) may be used for an early warning, thus complementing the geodetic ((Delta z)) monitoring and the (Lambda – t) plot (Fig. 11c).

Third, we consider the (h_{mathrm{w}}) and H dependence of the landslides. Our regime diagram (Fig. 3b) for (H simeq 20) cm showed that diverse slides occur when the (h_{mathrm{w}}) rises. The slide morphology includes 2 slides that progress downslope, spread, translation and flooding, which are not modeled by a standard stability analysis which assumes a slide with a simple geometry that occurs once. This result alerts that one should also consider the stability of the block that has already slid and the possibility of a basal slip. We found that the slide is triggered nearly simultaneous or after the groundwater seeps out from the cliff toe (Fig. 10c), an observation that may be used to constrain the lower limit of (T_{mathrm{slide}}). We showed that the (T_{mathrm{seep}}) is well explained using permeable flow velocity calculated using the permeability k formula in Eq. (11), indicating its applicability whenever direct measurement of k is not possible. We also showed that the (T_{mathrm{seep}}) increases with H (Fig. 10d), which we explained as a result of k decreasing with the effective stress (sigma ‘). This effect needs to be included when one seeks to improve the accuracy of estimating the seepage velocity and hence (T_{mathrm{slide}}). We also add that when the soil consists of fine particles, in order to improve the estimate of (F_{mathrm{s}}), one needs to include the capillary rise (h_{mathrm{c}}) and an uneven (sigma _{mathrm{z}}), latter of which is likely to be enhanced by cohesion.

Finally, we note the caveats of applying the experiments to field landslides. The first is the effect of friction, which scales in proportion to the (sigma _{mathrm{z}}) and, hence approximately to the H. This means that the small H in our experiments compared to field landslides disproportionately amplified the effect of cohesion (Iverson 2015), when the particle size consisting the cliff is the same. The experiments are not perfectly scaled models of field landslides. The second is the cliff geometry, whose height H to length L, and width W ratios were restricted in the range of (H / L = 0.30)–1.52 and (H / W = 0.11)–0.59. It is possible that this geometry, amplified the effect of arching and contributed to the (sigma _{mathrm{z}}) deviating from lithostacy. The third is the stepwise increase of (h_{mathrm{w}}) at (t sim 0) s, and the presence of an impermeable, smooth, rigid base. These conditions were suited to model seepage using a simple permeable flow model (Sect. 6.2), but nevertheless, it was simplified.