Temporal evolution of backward erosion piping in small-scale experiments

Backward erosion piping (BEP) is a form of internal erosion which can lead to failure of levees and dams. Most research focused on the critical head difference at which piping failure occurs. Two aspects have received less attention, namely (1) the temporal evolution of piping and (2) the local hydraulic conditions in the pipe and at the pipe tip. We present small-scale experiments with local pressure measurements in the pipe during equilibrium and pipe progression for different sands and degrees of hydraulic loading. The experiments confirm a positive relation between progression rate and grain size as well as the degree of hydraulic overloading. Furthermore, the analysis of local hydraulic conditions shows that the rate of BEP progression can be better explained by the bed shear stress and sediment transport in the pipe than by the seepage velocity at the pipe tip. The experiments show how different processes contribute to the piping process and these insights provide a first empirical basis for modeling pipe development using coupled seepage-sediment transport equations.

a Pipe depth A Pipe cross sectional area C u Uniformity coefficient, the ratio d 60 /d 10 D Depth of the set-up D* Dimensionless particle diameter D r Relative density (e maxe)/(e maxe min ) d x Particle diameter at which x% of sample (by weight) is finer e Void ratio = n/(1n) g Gravitational acceleration (

Introduction
Internal erosion is a significant threat to levees and dams [9,13,16,33], but time scales of the failure process are hard to quantify [12]. Different types of internal erosion can be distinguished: concentrated leak erosion, suffusion, contact erosion, and backward erosion piping. Backward erosion piping (BEP) can be defined as the failure process by which seepage under a structure erodes the granular foundation that is covered by a cohesive roof, forming a hydraulic shortcut. For many levees, backward erosion piping is the most common form of internal erosion because of the cohesive deposits which form a roof above sandy aquifers.

Piping process description
The BEP failure process starts with high external water levels that cause excess water pressure at the landward toe of the structure. When a cohesive blanket is present at the landward side, erosion can occur through either pre-existing defects in the blanket or cracks caused by uplift and rupture. If the vertical pressure gradient over the exit through the confining layer is sufficiently high to initiate vertical transport of sand, a cavity forms below the exit hole and the sand settles around the exit hole as a sand boil. At some point, a small pipe starts to develop from the cavity in the upstream direction. When the head difference (H) over the levee is sufficiently high and prolonged, the erosion continues in the upstream (backward) direction until a hydraulic shortcut forms, which likely results in a levee breach. However, if the head difference is too low or drops too quickly, an equilibrium occurs with a partially developed pipe.
The head drop at which equilibrium occurs (H eq ) varies with pipe length (l) [21,46]. The maximum H for which an equilibrium exists in stationary conditions is defined as the critical head difference H c , and l c is the corresponding pipe length. This point (H c, l c ) marks the transition from the regressive phase (l \ l c ) to the progressive phase (l [ l c ). Exceeding this head difference H c long enough ultimately leads to a levee breach.
Most previous BEP research focused on this critical head difference H c and the influence of aquifer geometry and sand properties on this critical condition, for example [4,21,46,52]. Such a stationary approach, which neglects the development over time, may be sufficient in many cases such as rivers with relatively long floods or dams with permanent pools. However, when the flood duration is short compared to the time scale of backward erosion, a dike may survive a short duration flood, whereas it would fail under a long duration flood [21]. These insights are important for levee safety assessments and for emergency response.
This raises the question of how to predict the temporal development of the backward erosion process. The twentieth century studies of Miesel, Müller-Kirchenbauer and Hanses [21,29,31] report very limited data on the development over time. More recently, several researchers studied pipe progression rates (dl/dt) experimentally [3,35,36,38,41,51]. Robbins et al. [38] used cylindrical flumes and correlated progression rates to seepage velocity and void ratio. Allan [3] studied progression rates under overloading (H [ H c ) in a rectangular setup with hole-type exit (seepage length L = 1.3 m). Vandenboer et al. [51] also studied the effect of overloading, but in a small-scale (L = 0.3 m) rectangular flume with hole-type exit. Pol et al. [35] derived progression rates from various available BEP experiments, and correlated these to hydraulic conductivity and the global horizontal hydraulic gradient. Robbins et al. [41] analyzed progression rates in small-scale flumes with a slope-type exit, for a wide range of grain sizes. Finally, Pol et al. [36] report pipe progression rates in a large-scale experiment. These studies [3,35,38,41,51] show that the progression rate is related to head difference, grain size and degree of compaction. Measured progression rates vary several orders of magnitude between different setups (e.g., slope or hole), sand types and degree of overloading.
A challenge in interpreting these results is that most experiments were performed on much smaller scales than typical levee dimensions (seepage length and aquifer depth in the order of 10-100 m). Since the progression rate appears to be a function of the applied and critical head, which do not scale linearly with seepage length, extrapolation to field conditions introduces a large uncertainty in predicted progression rates [35]. Scaling requires relating the progression rate to local, scale-independent conditions. Only Robbins et al. [38] and Pol et al. [36] evaluated local hydraulic gradients at the pipe tip, which are considered scale independent erosion criteria. The other studies only present global hydraulic gradients.

Piping process in non-equilibrium conditions
BEP consists of two distinct types of erosion [21]: lengthening of the pipe by detachment of grains from the soil skeleton at the pipe tip (primary erosion) and enlargement of the pipe cross section (secondary erosion). The primary erosion mechanism can be considered as successive slope failures which occur if the forces exerted on the grains by the seepage flow exceed the resistance [22,44]. Particles slide in the pipe, may rest temporarily on the bed, and are gradually transported by the flow. When detached grains rest temporarily on the bed or roll through the pipe, the shallower pipe results in a temporarily higher flow velocity, higher pipe gradient, and therefore a lower tip gradient (Fig. 1). That delays or stops the pipe lengthening. The tip gradient will increase due to the removal of grains in the pipe and the next slope failure occurs when the tip gradient has recovered to the critical value [37]. This intermittent pipe lengthening has been observed for example by Hanses [21]. The description above is based on observations when the pipe is close to equilibrium. In case of high progression rates, there will be continuously a layer of moving grains because the critical gradient at the pipe tip has been exceeded by such an amount that the influence of the sediment load in the pipe is not sufficient to reduce the pipe tip gradient below the critical value. In the case of internally unstable soils, the tip failures may also be delayed by erosion of the finer fraction from the coarser soil matrix [23]. In case of very fine dense sand, the rate of slope failures may be limited by dilatancy effects. Soil matrix expansion up to the critical porosity requires an inflow of water. The lower the permeability, the more time is needed to supply this water. Therefore it can pose a limit on the rate of primary erosion, similar to breaching flow slides described by Van Rhee [57]. At the same time, the pipe downstream deepens by secondary erosion if the bed shear stress exceeds its critical value [52]. Deepening reduces the gradient in the pipe and hence increases the tip gradient, so the two types of erosion are coupled.
Given these two erosion processes, the time scale of erosion can be considered as the combination of (1) time needed for erosion of grains in the pipe until the local tip gradient recovers to the critical value, (2) time needed for erosion of a finer fraction at the tip, and (3) time needed for dilation of soil matrix at the tip. For the sands considered in this paper, we assume that the second and third component can be neglected.

Modeling of pipe progression rates
Several methods have been proposed to model the development of the pipe length over time for engineering purposes, which are summarized here. Kézdi [24] hypothesized that the progression rate of BEP is proportional to the pore flow velocity at the pipe tip. However, his model of pipe progression neglects pipe resistance, seepage concentration at the pipe tip, and includes no critical pore velocity. Furthermore, there is hardly any experimental data to test this hypothesis. Some empirical relations between erosion rate and seepage gradient were developed in the context of streambank erosion of cohesive soils due to seepage [8,14]. Robbins et al. [41] determined such a relation for BEP experiments based on modeled seepage gradients and measured progression rates, and included this in a quasi-stationary BEP model. Additionally, several numerical models have been developed to predict BEP by coupling seepage, pipe flow and sediment transport relations [43,59]. However, the pipe length development as predicted by these models has not been validated on different types of BEP experiments, as experiments with suitable measurements are very limited.

Objective
The goal of this study is to better understand and model the development of piping erosion over time. To achieve this, we measured local, scale-independent, conditions that explain the temporal development of piping (i.e., the progression rate). We modified a commonly used BEP laboratory setup to measure pore pressures and pipe pressures during the piping process with a high spatial and temporal resolution. The experimental program included different sand types with varying degree of compaction to explore the effects of grain size and compaction. Part of the experimental method and a small part of the results were presented in Pol et al. [37]. The paper is structured as follows: chapter 2 describes the experimental method, chapter 3 the primary experimental results, chapter 4 an analysis of local flow conditions, pipe development and sediment transport, and chapters 5 and 6 contain the discussion and conclusions, respectively.

Modification of box-type setup
Previous research into BEP progression rates has either used a rectangular box-type setup [3,41,51] or a cylindertype setup [38]. The box-type setup poses less restrictions on flow from the sides and results in flow concentration toward the exit. On the other hand, a cylinder forces the pipe to grow right below a row of sensors to measure pipe pressures. To combine the advantages of both types, we modified the box-type setup used by [52] and [51] so that pipe pressures can be measured. The sample dimensions are 0.48 9 0.30 9 0.1 m (Fig. 2). The box has a 10 mm thick acrylate cover with a 6 mm diameter exit hole. The seepage length L equals 0.352 m, between upstream filter and the exit. The modification is twofold. First, two permeable barriers of filter fabric (0.05 mm aperture) were placed longitudinally and 35 mm apart to prevent sideward pipe growth but allow flow to pass. Second, silicon strips (0.3 mm high, 3 mm wide, 10 mm long) were placed diagonally with 20 mm spacing between the rows, and sand was sprinkled over them while the silicon dried (grey strips in Fig. 2). These two steps restrict the pipe path to the middle 35 mm of the box, without significantly influencing the flow. Hereby, the pipe can meander slightly, while also being close to the sensors and not significantly influencing the piping process. At the interface with the guides, the porosity may be slightly higher as the grains do not interlock. However, as will be demonstrated in Sect. 3.2, this does not have a significant effect on the results, as the pipe tip generally propagates at some distance from the guides. Finally, 20 mm spaced pressure ports were added in the center axis of the box and connected to pressure sensors at the side of the box. Influence of the ports on progression is expected to be negligible, given their limited volume and not extruding into the sand.

Materials and measurement techniques
The experiments include three fine, uniform sands with characteristics as shown in Table 1. The FPH sand was also used in Pol et al. [36]; the Baskarp B25 sand was used in Akrami et al. [2] and Rosenbrand et al. [42]. Grain sizes were determined by dry sieving. Porosity n min is based on the method in ASTM4253 (dry method, vibrating needle instead of vibrating table) and n max on ASTM4254 (funnel method). The porosity n sb and slope angle h sb of the sand boil per sand type were determined by measuring its height, diameter, and dry mass at the end of each test. The water temperature was between 20 and 22°C.
Water levels at the upstream and downstream sides were measured using riser tubes at least every 5 min, up to every minute during progression. The flow rate was measured at the same frequency using a digital scale. Pore pressures were measured by differential pressure transducers (Sensortechnics RPOP001D6A), through the ports P1-P19 (Fig. 2), at a sampling frequency of 10 Hz.
Three cameras recorded the erosion development. The main camera above the setup provides a top view of the sand sample, every 10 s. In some tests, we recorded short close-up videos of the erosion process using a second camera, which was placed temporarily on top of the cover. The last camera recorded the volume of the eroded sand, every minute. Pipe length and sand boil diameter were estimated visually at least every five minutes, up to every minute during progression (Fig. 3). Pipe depth in equilibrium conditions was measured using a laser device (DSE ODS 120) mounted on a movable frame to create transects perpendicular to the pipe. Pipe depth measurements during pipe progression were not feasible as the laser device would limit observation of the pipe by blocking the view from top and moving particles would affect the depth measurement.
Pipe flow velocities in equilibrium conditions were extracted from videos of dye tracer injections. A red dye was injected under low pressure through one of the pressure ports (usually P14 or P15) for approximately 2 s. The propagation velocity of the dye in the x-direction is obtained by tracking the change in color intensity relative to a video frame before injection [18]. Usually, we injected 5-7 times in the same equilibrium conditions, from which we take the median propagation velocity. We assume that this represents the maximum flow velocity over the velocity profile in a cross section (u max ).

Test procedure
The test procedure as described in [37] consists of sample preparation, loading, and measurement cycles. First, the sample is prepared with the box in the vertical position by sprinkling dry sand in de-aired water, and the sample is compacted by tapping the box with a hammer. Then, the box is closed and placed in the horizontal position, and the head at both sides of the sample is leveled. The loading procedure is as follows: keep the head difference constant if there is still erosion after 5 min, or increase it otherwise. The head difference is increased by reducing the downstream head, and the upstream head is kept constant. When the pipe reached the upstream filter, the downstream head was raised (suddenly) to stop erosion. Then, the head difference was increased in small steps of 2 mm until grains started moving somewhere in the pipe. If the movement continued, the head difference was decreased by 2 mm. This procedure was iterated until the bed was just in equilibrium. In those conditions, the head drop was kept constant and the pipe depth, local hydraulic gradients in the pipe, and local maximum flow velocity u max were measured to determine the critical bed shear stress. During the test, measurements of pipe length, sand boil radius, flow rate and total head drop were taken every five minutes, up to every minute during progression. After the test, the sand boil was collected, dried and weighed.

Test program
The test program consists of two phases as shown in Table 2. First, tests 217-222 are reference tests without guides to verify that the changes in experimental setup do not influence the critical head difference and average progression rate. Note that tests 218 and 219 are not representative as the pipe initially developed toward the side of the box, affecting the critical head. Tests 223-231 were not included in Table 2 since these involved iterative improvements of the guides to obtain stable results. Phase 2 consists of tests 232-250 with the adapted setup and varying sand type, degree of compaction and hydraulic loading. Densely packed samples have a relative density (D r ) of 0.7-0.8 and loose samples of 0.5-0.55. The standard loading scenario L1 is to gradually increase the head difference to the critical head difference (H c ) and then keep it constant. L3 and L4 are overloading scenarios, which means that after a stable pipe developed with pipe length l & 0.14 m (Table 2), the head is suddenly increased to 1.2ÁH c (L3) and 1.1ÁH c (L4), respectively, and then kept constant. Note that the effective head drop over the sample (H corr ) is not always constant after H c is reached, due to filter and exit losses changing with flow rate. Correction for filter and exit losses is described in Sect. 3.2. Furthermore, the hydraulic conductivity k is estimated from flow rate and pressure gradient near the upstream filter, so it may be less reliable for tests 217-222 (marked with #). The results in Table 2 are discussed in Sect. 3.2.

Experimental results
This section reports the experimental observations and the basic measurements of hydraulic head and pipe geometry. Analyses of progression rates, shear stresses and sediment transport require more interpretation and are therefore described separately in Sect. 4.

General observations on the erosion process
This section describes the observed phases of the BEP process and (visual) observations of the erosion process at the grain scale. Like in other hole-type experiments [29,51,52], each test showed several phases for an increasing head drop: fluidization of sand in the exit hole, formation of a circular void (lens) around the exit hole, pipe growth toward an equilibrium (regressive) and progressive pipe growth until the pipe forms a hydraulic shortcut. Figure 4 indicates these phases in a plot of the head at each transducer. Before the erosion phases, the heads respond almost instantaneously to an increasing head drop. During erosion, the heads decrease also gradually under a constant head drop due to the pipe lengthening. It was also observed that after a head increase, the location in the pipe where grains started eroding was not always the same. Visual observations during the head increase indicate that erosion sometimes starts in the bed, sometimes at the tip. This suggests that both the bed and the tip are close to critical conditions, at least in the regressive phase. At a few instances during the test program, the erosion process at the pipe tip was visualized using close-up videos. Grain detachment at the tip generally occurred in cycles. Sometimes, first a small displacement of particles in the zone upstream of the tip was observed, which increased the porosity locally. Sometimes, there was rearrangement of a few small particles in the sand upstream of the tip without the other grains moving. Quickly after the small displacement, a group of grains detached and moved into the pipe. Part of the grains washed away directly, while another part settled close to the tip. These settled grains were transported gradually until the cycle repeated. Between the group detachments, also individual particles detach. At higher progression rates, it was more difficult to distinguish separate cycles and the erosion process is more continuous. There is much variation between and within tests regarding the occurrence and duration of these steps in the erosion cycle, but the process generally followed this cycle.
One test (B25_245) includes a close-up video of the tip during a transition from equilibrium to erosion and simultaneously the pipe tip grows closely under the pressure sensors, which allows to observe the effect of particle detachment on the pressure gradients at the tip and in the pipe. Based on this example, Pol et al. [37] show that the pipe resistance induced by the detached particles reduces the tip gradient temporarily below the critical tip gradient and temporarily stops the tip erosion. These observations indicate that the transport in the pipe affects the progression rate.

Critical head, pipe length and hydraulic gradients
The main results of the experiments are given in Table 2.
The measured critical head drop H c is the sum of the head drop over the sample H c,corr , the upstream filter loss, and the exit losses. The filter loss was estimated by a linear regression of the head profile through sensors P14-P15. The exit loss is estimated from regression on sensors P2-P5, but due to nonlinearities in the head profile around the tip this is only reliable when the tip has passed P5 (x [ 0.25 m, l [ 0.115 m). Therefore, for l [ 0.115 m, the exit loss from regression was related to the flow rate Q for each test. Ultimately, the estimated exit loss throughout the test is based on this relation between exit loss and measured flow rate. At the critical pipe length, the average estimated filter and exit losses were 6.1 and 5.6 mm, respectively. The critical pipe length l c = x tip,c -x exit is the pipe length when the head reaches H c . The average progression rate after the critical head has been reached v c,avg = (L-l c )/(t end -t c ) is based on visual observations of the pipe tip position. Hydraulic gradients between transducer pairs were derived from the pressure measurements, mostly having a 0.02 m spacing (Fig. 5). Critical tip gradients i c,tip are defined as the maximum gradient of the transducer pairs passed during the regressive phase (so pipe is close to equilibrium). These values are not available for all tests, as the pipe tip sometimes passed besides all pressure ports in the regressive phase. If multiple values were obtained for one test, Table 2 gives the maximum. Modification of the setup (guides and pressure ports) did not lead to significant differences between the control group of dense B25 sand (tests 217 and 220-222, see Table 2) and modified group (tests 232, 236, 243), in terms of critical head H c,corr (6.3-6.8 cm before, 5.9-6.2 cm after), critical length l c (10-19.5 cm before, 15.6-18.6 cm after) and progression rate v c,avg (0.05-0.09 mm/s before, 0.05-0.10 mm/s after). The slightly lower H c (9%) after modification may be the result of a locally lower D r and hence higher conductivity at the guides. In case of the coarse FS35 sand, however, the pipe width may have been limited by the guide distance of 35 mm (see Sect. 2.1), resulting in slightly deeper pipes compared to a situation without guides.
Both the net head drop over the sample H c,corr (Table 2) and the local critical tip gradient i c,tip have a positive relation with D r for both the B25 and FS35 sands (Fig. 6), which confirms findings by [38]. The critical pipe length l c varied between 30 and 60% of the seepage length. Given the large variability in l c , it is not possible to infer significant effects of grain size or compaction.

Pipe length development
Kezdi [24] expected acceleration of the pipe development because of the increasing upstream secant (average) gradient with increasing pipe length. Figure 7 shows the pipe length development between the critical length l = l c and when the upstream filter is reached (l = L). Pipe length and time are normalized. The normalized pipe length becomes , in which t c is the time when l = l c and t end is the time when l = L. On average, there is some acceleration in normal loading tests, which are closer to equilibrium. However, there is hardly acceleration in the B25 overloading tests as indicated by the nearly linear curve (note that absolute progression rates are higher).
This confirms experiments by Vandenboer et al. [51] using the same setup, and another hole-type experiment by Miesel [30]. However, slope-type experiments by Robbins et al. [41] showed very rapid progression and did not accelerate. Slope-type tests are initiation-dominated [52] and therefore overloaded more severely under a constant head. Apparently, the progression rate does not increase with pipe length in those overloading conditions, despite the increasing upstream secant gradient. That indicates that the progression rate is limited by the transport of sediment down the pipe rather than the limited supply of water to the pipe tip in severely overloaded conditions. The observation of acceleration shows that the progression rate (load effect) is not constant for a constant load. Here it is noted that the upstream filter resistance increased slightly with pipe length in several tests, resulting in a decreasing head drop over the sample. Without that resistance, the acceleration is expected to be even more pronounced.

Pipe geometry after test
The pipe geometry was analyzed to estimate the shear stress acting on the pipe bottom during equilibrium in a fully formed pipe. When a pipe had fully developed to the   When the pipes progressed, the channels were often migrating sideward, especially in the downstream parts of the finer sands B25 and FPH. This results in inactive channels without sediment transport and presumably a low flow rate. As we assume that the main flow conveying channel is representative for the pipe flow conditions such as shear stress, using the entire channel would lead to an incorrect bed shear stress s (as a avg is underestimated). This channel migration is not the result of the setup modification as it also occurs in [50], though it is expected to occur less due to the application of the guides, see Fig. 2. The definition of the boundaries of the main channel is based on visual interpretation of the depth profile (Fig. 8). If one channel is clearly larger than the other, the largest is selected as main channel. If both are equally large or there is only one, the main channel equals the entire channel. Figure 9 shows pipe geometry based on the entire channel and the main channel. The main channel approach yields higher a avg , lower w avg and more consistent geometries across the three sands. In Sect. 4.1 (Fig. 10) we show that the main channel approach yields more consistent shear stresses. Therefore, the main channel geometry is a better representation of the pipe flow than the entire cross section and we use it to estimate pipe flow conditions. Figure 9 shows the resulting depth profile (a avg ) per sand type. The solid line indicates the average of tests 232-250, the error bars indicate ± one standard deviation. The pipe depth scales with the permeability as a $ ffiffiffi j 3 p [50]. Scaling with d 50 gives reasonable results too for these sands; which one scales better can only be assessed with a wider range of grain size and uniformity. The depth profile can be approximated with a power function as: The profiles of a max and R (not shown in Fig. 9) have a similar form with coefficients 10 and 2.8, respectively. The shape of the depth profile will depend on the spatial distribution of seepage toward the pipe, so on aquifer geometry. As this relation is purely empirical, it is only to be used to analyze the current experiments. Three B25 tests (230, 236 and 245) include depth measurements when the pipe was partially developed (l \ L). From the depth data (not shown here), it follows that the partially developed depth profile is similar to Eq. 1. Therefore, we assume that Eq. 1 also holds during progression. Figure 9b indicates that there was no significant difference in equilibrium depth profiles between the average of all tests and the ones with overloading (L3 and L4). Therefore, we assume that overloading does not lead to significant differences in depth.
The average pipe width w avg is rather constant along the pipe (Fig. 9d) and in the order of 60 times d 50 for B25 and FPH sands. This results in a main channel w avg /a avg ratio of 10 at the downstream side to 25 near the pipe tip. Furthermore, based on the photos taken from top, we found the tip width to be approximately 30 times d 50 , confirming [53,63]. However, the tip width increases with the degree Fig. 10 Critical Shields number during equilibrium with fully developed pipe in tests 232-250, average value per sand type, plotted with empirical relations [6,52] and experiments in laminar flow in classical flumes [19,25,28,34,55,58,60,61,65] and piping cylinders [55] Acta Geotechnica (2022) 17:4555-4576 4565 of overloading, up to 45Ád 50 at 1.2ÁH c . The pipe width in FS35 sand is lower than in B25 and FPH sand, with w avg = 40Ád 50 . This may be caused by the limited space (35 mm) between the filter fabric guides. This restriction may have caused narrower and deeper cross sections, which could explain the underestimated FS35 depth using Eq. 1 (Fig. 9).

Pipe progression analysis
This section analyzes the pipe progression rates and its relationship with tip seepage velocity and pipe shear stress. First, the bed shear stress in equilibrium conditions with a fully developed pipe is computed with different approaches and compared to existing data of critical bed shear stress. This information is subsequently used to determine the bed shear stress during pipe progression. The tip seepage velocity is calculated from the measured tip gradient i tip , conductivity k and porosity n, and both parameters are related to the pipe progression rate. Finally, the sediment transport is compared to sediment transport relations from classical, laminar flume experiments conducted for studies of sediment transport in open channels.

Critical bed shear stress measurements during equilibrium
The critical bed shear stress for incipient motion is an important parameter to predict piping and was determined from measurements at the end of each test, for the situation of a fully developed pipe. The pipe geometry is based on the main channel (see Sect. 3.4). The local hydraulic gradients between transducer pairs are only calculated for those pairs below which the pipe progressed. The maximum flow velocity in the pipe was measured by injection of a dye in one of the pressure ports close to the filter (P14 or P15). Bed shear stress can be calculated based on either pressure gradient and pipe depth, pressure gradient and flow velocity or flow velocity and pipe depth. First, we present the three equations to calculate shear stress. Later, we show the resulting shear stresses for the piping experiments. In the equations below, the subscript of s indicates on which quantities it is based, for example s ai is based on pipe depth and hydraulic gradient. First, from a balance of forces it follows that the average shear stress along the wetted perimeter is the product of hydraulic radius R = a/2 and local pressure gradient in the pipe dp/dx, which gives for parallel plates [46]: Alternatively, one can use the relation between depthaveraged flow velocity U, pressure gradient and pipe depth from the Poiseuille equation for laminar flow between parallel plates [46]: dp dx In which U ¼ 2=3u max [48]. Substitution of a in the shear stress Eqs. (2) by (3) gives: Alternatively, substitution of dp/dx in the shear stress Eqs. (2) by (3) gives: Note that the Reynolds number (Re = 4RU/v) during equilibrium at the end of the tests was in the order of 40 for FPH sand, 80 for B25 sand and 280 for FS35 sand, which confirms that the flow was laminar (Re \ 2100). As the main channel w/a ratio lies in the range of 10 to 25, the assumption of parallel plates is a reasonable approximation [48,52].
These three expressions (Eqs. 2, 4 and 5) for the shear stress are equivalent under the assumption of laminar flow between parallel plates but require different input. During equilibrium conditions, all three inputs (pipe depth a avg , local gradient between transducer pairs i, and flow velocity u max ) were measured along the pipe to be able to compare Eqs. 2, 4 and 5. Finally, the measured shear stress s(x) is averaged over the length of the pipe. In the analysis we assume that water density q w = 1000 kgÁm -3 , and dynamic viscosity l = 1 mPaÁs (20 degrees C). Figure 10 shows the resulting critical bed shear stress from Eqs. 2

Drivers of the progression rate
The progression rate v = dl/dt is a practical engineering metric for the temporal development of the piping process. This section correlates the observed progression rate during the progressive erosion phase to two variables: (1) seepage velocities upstream of the pipe tip following Kézdi [24], and (2) bed shear stresses in the pipe. The seepage velocity just upstream of the pipe tip is derived from measured tip gradients: In which i tip is the local hydraulic gradient measured over a distance of 2 cm (transducer spacing) upstream of the pipe tip, and k and n are the initial hydraulic conductivity and porosity. While Kézdi [24] neglects head loss in the pipe and assumed that i tip equals the average upstream gradient, this section applies the same concept, except for using measured local gradients. In contrast, [41] used backcalculated tip gradients from a FEM model.
The bed shear stress is calculated with Eq. 2 during the progressive erosion phase, since no flow velocity measurements during progression are available: sðx; tÞ ¼ 0:5q w gaðx; tÞiðx; tÞ ð 7Þ The pipe gradient i is calculated as the hydraulic gradient between the transducers downstream of the pipe tip. Shear stress is calculated both based on the average of all transducer pairs downstream of the tip (s bed;average;all ) and based on the average of the transducer pairs where the pipe passed right under both ports (s bed;average;passed ). The first may include more scatter from transducers pairs which are partly above and partly besides the pipe but contain more datapoints to be included in the averaging.
Due to a lack of depth measurements during progression, we assume that the depth profile has the same shape as during equilibrium at the end of the tests (Eq. 1): The depth of 0.0001 m is subtracted to account for the residual erosion after reaching the upstream filter (see Sect. 3.4). We believe this provides a reasonable estimate of the depth profile as it is similar to some depth measurements made for partially developed pipes in three B25 tests (see Sect. 3.4). Furthermore, Eq. 8 was verified for test B25_247 (at l = 0.14 m) using the Poiseuille relation (Eq. 3) and measured gradient near the exit (i = 0.153 between transducers P2-P4), flow rate (Q ¼ Uwa = 1.12 mL/s), and pipe width (w = 0.01 m). Equation 3 yields a calculated depth of a = 0.96 mm, whereas Eq. 8 yields a = 1.1 mm.
We calculated seepage velocity (Eq. 6) and pipe shear stress (Eq. 7) for each time that the tip passed a pressure port in the progressive phase (l [ l c ). However, we omit seepage velocity and tip gradient data if the pipe tip passed besides that transducer, as this data are unreliable. And we omit shear stress data if no value is measured for either s bed;average;passed or s bed;average;all . This results in 75 datapoints from test 232-250 for seepage velocity and 103 datapoints from test 233-250 for shear stresses, which are from different sand types and degrees of overloading and measured at different pipe lengths. Note that using an equal number of datapoints for all parameters in Fig. 11 would not affect the trends in Fig. 11, but the lesser datapoints would make the analysis in Table 3 unreliable. Finally, the corresponding progression rate at these passing moments is calculated using a moving average (over 3 datapoints; usually 3 min) of the visually observed tip position. Figure 11 shows the relations between progression rate and several parameters related to seepage velocity and pipe shear stress. Values of q p indicate Pearson correlations. Figure 11a shows that the progression rate is proportional to the tip seepage velocity u p for tests with normal loading (L1) on different sand types. However, the overloading tests (L3 and L4) on B25 sand show higher progression rates which cannot be explained by the seepage velocity as the seepage velocity hardly changes under overloading. Figure 11a-c show how the effect of u p is composed of k and i tip and that the hydraulic conductivity explains most of the variation in progression rate observed in normal loading tests, not the tip gradient. Furthermore, the measured local tip gradient i tip does not seem to increase with the degree of overloading, although there is substantial variability in these measurements. This can be caused by a higher porosity at the tip, or by higher resistance in the pipe (due to high sediment load). Furthermore, it is likely that the tip gradient cannot exceed its critical value much, as the tip material will collapse (see Fig. 1). Figure 11e-f shows that the bed shear stress is a better predictor, see also Sect. 4.3 for more discussion. Compared to the seepage velocity it predicts the overloading tests reasonably well and explains part of the variation within groups (e.g., group FS35,L1). s bed;average;all gives an even stronger relation (Fig. 11f), probably due to the averaging of noise caused by some pressure ports being besides the pipe and that generally a larger part of the pipe is included in the averaging due to less strict criteria. Finally, the data show no relation between D r and progression rate in normal loading tests (Fig. 11d). We fitted a linear regression line through the same data as in Fig. 11e and f, of which the intercept is a proxy of the critical bed shear stress during progression. The same procedure was followed for the shear stress at the first transducer pair downstream of the tip and for the average shear stress between the five transducers (P2-P6) closest to the exit. Table 3 shows the obtained results, including the values with fully developed pipes at the end of the tests. It appears that the estimated s c during progression is lower than during equilibrium in fully developed pipes. The estimated s c near the exit is similar to the average s c in the pipe, but the estimated s c near the pipe tip is clearly lower. There is a general consensus that upward seepage reduces both the occurring shear stress due to a change in velocity profile as well as the critical shear stress [7,15,26]. Regarding the impact on sediment transport rates, these studies show diverging results. Upward seepage effects were omitted from the Sellmeijer model [47] as they expected that this would not affect the eroding grains on top of the bed. As upward seepage is most severe at the pipe tip, this effect could explain the lower critical shear stress near the tip in Table 3. The critical bed shear stress under upward seepage is equal to s c ¼ s c;0 1 À i u i f where s c,0 is the critical bed shear stress without seepage, i u is the upward hydraulic gradient in the bed, and i f is the critical gradient for fluidization [7]. Finite element simulations of test B25_245 (for l = 0.185 m) [39] show that the upward seepage gradient below the pipe varies from approximately 0.35 at the tip to 0.2 at a distance of 3 cm from the tip. The fluidization gradient is i f ¼ ð1 À nÞ q s =q w À 1 ð Þ [49], which equals 1 for B25 sand with n = 0.4. The resulting 20-35% reduction in s c partly explains the lower experimentally obtained critical shear stress 'during progression' as compared to 'during equilibrium' (Table 3). These findings indicate that upward seepage may have a significant effect on s c in BEP. Upward seepage affects the acting shear stress s through a change in the near-bed velocity profile [7]. A full quantification of these effects for the case of laminar flow in rectangular ducts is beyond the scope of this research and needs to be verified with modeling or future experiments with more detailed information on pipe depth, pipe velocity profile, and upward seepage gradients in case of partially developed pipes. Omitting the upward seepage effect from BEP models may be acceptable if the effect on the acting shear stress s balances the effect on the critical shear stress s c .

Sediment transport
The dependence on the bed shear stress raises the question whether the piping erosion rate can be predicted using sediment transport equations. Sediment transport (bed load) in turbulent flow is typically related to the bed shear stress through: the coefficients for turbulent flow are a 1 % 4; a 2 % 1:5 [62]. Similar equations exist for laminar flow [5,6,27,32,64]. The sediment transport rate in laminar flow scales better by a viscous scaling rather than the inertial scaling of the Einstein number [5,32]: Although our experiments were not specifically designed for this purpose, we complemented the available experimental data on sediment transport rates in laminar flow [1,5,10,11,20,27,45] with estimated sediment transport rates from the piping experiments in this paper. Table 4 summarizes the materials and flow conditions in the classical flume experiments. The sediment transport rate q v (volume per unit width) in the piping experiments is Based on linear regression on the data of progression rate versus shear stress (Fig. 10). # = not significant (p [ 0.05): s crit estimated visually. * = less significant (0.01 \ p \ 0.05) Acta Geotechnica (2022) 17:4555-4576 4569 estimated from the growth rate of the sand boil radius (Fig. 3) and an average main channel width near the exit of 40*d 50 (Fig. 10). Assuming a cone-shaped sand boil with equal inner and outer slopes, its sediment volume is given by: In which n sb the porosity in the sand boil, h sb the slope angle, r b the radius of the sand boil (observed during the tests), and r e the exit hole radius (0.003 m). Values of n sb and h sb given in Table 1 are based on sand boil dimensions and dry sand mass measured at the end of tests 224-250. The Shields number H at the downstream end of the pipe is based on Eq. 7 using the average pipe gradient over the 5 transducers near the exit (P2-P6) and pipe depth a avg at the exit (x = 0.125 m) from Eq. 8, which increases with pipe length. The critical Shields numbers H c at the downstream end of the pipe are calculated as 0.0736 (FPH), 0.0705 (B25) and 0.0615 (FS35) from s exit in Table 3, so from extrapolating the shear stress to v = 0. Note that this estimate of H c during progression is relatively low compared to the value at the end of the tests (Fig. 10), which is possibly due to upward seepage during the piping experiments. Figure 12 shows the sediment transport rate q Ãv as function of H H À H c ð Þ , for both the classical flume experiments and the piping experiments. The experimental data mostly follows the same trend. Only the data of Charru et al. [5] is below the range of other data, and Grass and Ayoub [20] shows a higher exponent. The transport rate in the piping experiments is approximately a factor 5 below the average trend, but falls within the experimental range, and shows a similar trend with increasing Shields number. As can be expected, the scatter in the piping experiments is larger than in classical flume experiments because there is intermittent erosion, no wide uniform bed and flow profile, and no depth profile measurements during progression. Fig. 13 Experimental data in classical flume experiments (see Fig. 12) as predicted by several sediment transport equations including Eq. 12 Least square fitting of the flume data alone results in the following empirical relation with R 2 = 0.90 and a median absolute percentage error (|q predictq exp |/q exp ) of 0.373, drawn in Fig. 12: Figure 13 shows the ability of Eq. 12 and the models by Yalin [64], Charru et al. [5], Cheng [6] and Ouriemi et al. [32] to predict the classical flume experiments from the literature [1,5,10,11,20,27,45]. Cheng [6] and Ouriemi et al. [32] fit well to the higher transport rates, but overestimate the transport for flow conditions close to critical. Yalin [64] generally performs well but shows some underestimation closer to critical. The empirical model by Charru is about one order of magnitude below the average trend of the flume data, which is in line with findings of Wewer et al. [59] that Charru underestimates the pipe progression. However, it is a lower bound for the piping data. The simple empirical regression from Eq. 12 performs well over a wide range of measurements but underestimates the highest transport rates. As in the piping experiments q Ãv \ 10 -2 , Eq. 12 is a suitable empirical model.
So far, the analysis focused on the relation between transport rate and shear stress in the classical flume experiments and piping experiments. The measured transport rate can also be used for an additional cross-check on the assumptions on the pipe geometry during progression. For a progressing pipe, the sediment mass balance states that the increase in pipe volume per unit time (dl(1n)ÁA exit ) is equal to the total sediment transport at the downstream end of that pipe (Q v Ádt). Assuming a rectangular cross section (A exit = w exit Áa exit ), this relates the progression rate dl/dt to the sediment transport rate as: Combining Eqs. 12 and 13 results in a progression rate as function of the measured shear stress near the exit: For the same data as in Fig. 10 and using a exit from Eq. 8 and assuming w exit = 40Ád 50 , Eqs. 13 and 14 were evaluated and compared to the measured progression rate. Figure 14a shows that v predict;V yields a correct magnitude of dl/dt, providing a cross-check that reasonable assumptions were made regarding a exit and w exit . Figure 14b shows that Eq. 14 results in an over-prediction of the progression rate by a factor 5, which agrees with the difference between piping tests and classical flume experiments in Fig. 12. Therefore, we recommend using a coefficient of 0.08 instead of 0.39 in Eq. 12 when modeling BEP using sediment transport relations, until more is known about the applicability of these sediment transport relations for BEP.

Discussion
This discussion section compares the results to previous experiments, reflects on the experimental setup and finally discusses the erosion process and some implications for modeling the temporal development of BEP.

Comparison with other experiments
First, we compare our results with previous BEP experiments. No experiments with progression rates on these specific sands and with the same setup have been reported, but there are several experiments with a similar setup and sand of a similar size. The critical head of B25 is similar to tests on M32 sand (d 50 = 0.251 mm) by Vandenboer [50], but is lower than small-scale scale tests on Itterbeck due to lower permeability. We find a slightly larger pipe depth at the end of the tests compared to Vandenboer [50], but this can be explained by differences in test procedure, as Vandenboer stopped the experiment immediately when reaching the filter, whereas we brought the pipe in equilibrium conditions. In terms of progression rates, the closest comparison can be made between the tests on FPH sand and tests on M34 sand (d 50 = 0.155 mm) by Vandenboer et al. [51], which were both conducted in small-scale setups with hole-type exit. For loading conditions of H/H c = 1 (like L1), Vandenboer reports an average progression rate of 7.3Á10 -5 m/ s which is comparable to our rate of 4.7Á10 -5 m/s for FPH sand. Comparing with small-scale slope exit experiments by Van Beek [35,54] and Robbins et al. [41], we find almost one order of magnitude lower progression rates. As explained in Sect. 3.3, this is likely due to fact that slope setups are generally initiation-dominated and therefore are more severely overloaded when the pipe length increases. This indicates that a model with local parameters can be more accurate compared to one using average gradients and average progression rates. Robbins et al. [41] proposed a relation with the local tip gradient: In which e = n/(1n) is the void ratio. Figure 15 shows that Eq. 15 cannot predict our measurements accurately for both normal loading and overloading based on measured gradients.

Reflection on the analysis and experimental setup
In this section, we discuss several assumptions in our analysis which may influence the results. First, this research focuses on the internal erosion mechanism of backward erosion piping, which can occur when there is a supporting cohesive roof above the cohesionless soil. In this case, erosion takes place predominantly in the sand, resulting in relatively wide pipes which can be approximated by flow in wide rectangular ducts. In case of other forms of internal erosion, such as concentrated leak erosion, the erodible surface is all around, resulting in circular pipes. Findings in this paper cannot be applied directly to other forms of internal erosion. Due to the assumption of laminar flow, the sediment transport equations may also not apply to the process of pipe widening (enlargement) after a hydraulic shortcut formed, which likely results in turbulent pipe flow.
Second, the critical local tip gradient (Fig. 5b) can serve as a local criterion for pipe progression. Like a cylindrical setup [38], our adapted setup is very suitable to measure these criteria. However, it is important to note that its value depends on the sensor spacing, and when applied in numerical modeling this should match with the grid size for instance.
Third, as we have no frequent pipe depth measurements during progression, we assume that the shape of the equilibrium depth profile at the end of the tests also applies to partially developed pipes. This assumption affects the computed bed shear stresses. Yet, we expect this provides a reasonable estimate of the depth profile as it matches depth measurements made for partially developed pipes in three B25 tests 230, 236 and 245.
Fourth, based on a preliminary analysis we found a difference between the critical bed shear stress in equilibrium conditions at the end of the tests and during progression. Furthermore, during progression, the estimated critical shear stress near the pipe tip is lower than near the exit. A possible explanation for this finding is the effect of upward seepage on the grain stability, which is absent the end of the tests. Upward seepage also affects the pipe flow velocity profile, and thus the relation between pipe flow rate and bed shear stress. For modeling purposes, it is important to apply a critical shear stress that is consistent with the applied pipe flow equations. It is recommended to investigate whether computations with and without upward seepage give a significant difference in critical pipe gradient.
Finally, we find that the sediment transport rate in our experiments is approximately a factor five below what is expected based on most classical flume experiments, Fig. 15 Measured progression rate in piping experiments and prediction using Eq. 15 [41] although still above Charru's data [5]. The reason for the lower transport rate in piping tests is not clear, as there are several possible explanations. First, the ratio of the flow depth to the particle diameter (a avg /d 50 ) is relatively low, which may restrict the transport. Second, the intermittent erosion in BEP may result in a lower transport rate as the change in sand boil size is measured over a longer interval (typically 5 min). Furthermore, values of H and H c during progression depend on the assumed depth profile, but this cannot explain a factor 5. Finally, the grains may be transported over a smaller width than 40*d 50, concentrated in the channel center, which would increase q Ãv , but also cannot explain a factor 5. Given the uncertainties in the determination of H, H c and q Ãv during progression, we recommend testing whether this difference can be confirmed by BEP experiments with more detailed measurements of sediment transport rate and both acting and critical shear stress.
The adapted setup allows to force the pipe development in the vicinity of the pore pressure sensors, and thus allows to measure local pipe gradients and flow velocities. This in turn, allows to determine local hydraulic gradients and bed shear stresses during equilibrium and pipe development. We see some further improvements of the setup to obtain better measurements. First, additional pressure sensors to measure the exit head loss and vertical gradients below the pipe. Pressure ports using wider slots instead of small holes would allow for a more accurate measurement of tip gradients, as it reduces the scatter from the pipe passing at some distance from the center line. Second, to reduce uncertainties in the pipe depth development during progression, more geometry measurements should be taken during the regressive and progressive erosion phase. This will yield more reliable estimates of acting and critical shear stress during pipe development.

Erosion process
Pipe lengthening occurs if the stability criterion at the pipe tip (e.g., local gradient) is exceeded [21]. On the other hand this often requires pipe deepening which is governed by the critical bed shear stress in the pipe [46]. As conceptualized in the introduction, the time scale of the BEP erosion process (time between successive tip failures) can be dominated by gradual removal of finer material from the soil matrix at the tip, by gradual soil matrix expansion due to limited inflow (dilatancy), or by the sediment transport capacity of the pipe. In these uniform sands and in nearcritical conditions, our experiments show no gradual removal of fines from the tip, as assumed by Fujisawa et al. [17] amongst others, but instead a stable soil matrix at the tip between successive tip failures. If the process would be driven by gradual removal of fines or by dilatancy, one would expect a strong relation between progression rate and tip gradient, which was not found in both normal and overloading conditions. It was also observed in near-critical conditions that grain detachment leads to a temporary increase in pipe gradient and corresponding decrease in tip gradient, indicating that the time between tip failures depends on the time needed for clearing the pipe. This importance of the sediment transport capacity is supported by the strong correlation between progression rate and bed shear stress in our experiments, in both normal and overloading conditions. We expect that in practice the pipe sediment transport capacity dominates the BEP time scale, because most river levees susceptible to piping are built on uniform sand and are only moderately overloaded in terms of critical head. Removal of finer particles may be important for graded material, and dilatancy may start to play a role in case of very strong overloading of dense sand.

Implications for modeling
The results indicate that the progression rate is more complex than just a function of the seepage velocity or gradient at the tip. During overloading, for example, the progression rate increases much more than is expected based on the seepage velocity only. Because the sediment transport rate in the piping experiments shows a similar relation with the bed shear stress as in classical flume experiments in laminar flow, we propose to use a transport relation like Eq. 12 for BEP modeling, similar to the recent work by Wewer et al. [59]. However, as progression occurs on exceedance of a critical condition at the pipe tip [40], such a condition would further improve the model. Therefore we propose to implement a sediment balance (Exner), a transport relation like Eq. 12 and a critical tip gradient in a 3-dimensional BEP model [56]. Such a coupled model may still need calibration on experiments like the ones presented here. Finally, we note that values assumed for H c must be consistent with the modeling of H, in terms of whether is accounted for effects of vertical seepage on the velocity profile and the bed shear stress.

Conclusions
To better understand the processes determining the temporal development of backward erosion piping, we performed a series of small-scale experiments. Our experimental setup guides the eroding pipe along a densely spaced row of pressure transducers, which allows local pressure measurements right in the pipe during equilibrium and progression. This is an advantage compared to most previous experiments [3,41,51]. Pipe geometry and flow velocity measurements allow for calculating bed shear stresses. The guiding of the pipe did not significantly influence the progression rate. We used three fine uniform sands (0.185 \ d 50 \ 0.422 mm), compacted to a relative density between approximately 50% and 80%. In addition to the regular loading up to the critical head, some experiments were overloaded at 10-20% above the critical head.
The results confirm the positive relation between progression rate and grain size or permeability [3,35,38], as well as the distinct effect of even a small degree of overloading [51]. The finding that the progression rate depends on the loading introduces an important problem. For a given head drop H, the loading effect which drives progression (e.g., tip gradient or bed shear stress) depends on both the scale of the problem (i.e., seepage length) and the length of the pipe. Therefore, the progression rate will be time-and scale-dependent. Kézdi [24] and Robbins et al. [41] proposed to relate the progression rate to the local hydraulic gradient or seepage velocity at the pipe tip, which is modeled as function of time and scale. However, the current experiments indicate that the sediment transport in the pipe is a stronger predictor than the seepage velocity at the tip. First, the tip gradient does not increase significantly with the degree of overloading. Second, groups of detached sediment create extra pipe flow resistance and temporarily reduce the tip gradient, and thus control the time scale of tip erosion. Third, the estimated bed shear stress during pipe progression correlates well (q = 0.84) with the progression rate, compared to q = 0.33 for seepage velocity. Especially the overloading situations are predicted better by bed shear stress than by seepage velocity. This suggests that the sediment transport capacity is a main, limiting factor in BEP progression rates of uniform sands.
The critical bed shear stress has been determined for the first time in a rectangular piping setup. The values after the tests (complete pipe) are in good agreement with piping tests in a cylindrical setup [55] and classical flume tests in laminar flows. However, the estimated critical shear stress for a partially developed pipe are much lower and seem to decrease in the vicinity of the pipe tip. This reduction can be partly explained by the effect of strong upward seepage [63]. Seepage is known to decrease both the exerted and critical shear stress in turbulent flows [7], but this has not been fully quantified for laminar pipe flow, and requires further study.
The relation between sediment transport rate and Shields number in the BEP experiments follows a similar trend as in classical, laminar flume experiments [e.g., 20]. We fitted an empirical relation (Eq. 12) to describe the sediment transport as function of the Shields number. Since there is a large scatter in the data from both the piping experiments and flume experiments at low transport rates, this should be studied further with more accurate measurements. If these relations indeed apply to backward erosion piping, timedependent pipe progression can be modeled by coupling a BEP model [56] with sediment transport relations (Eq. 12) and a sediment mass balance [59] and a local critical tip gradient as primary erosion criterion [40]. Such a model would contribute to more accurate assessments of the risk of piping failure in case of a limited duration of the hydraulic load, and therefore support more efficient designs of flood defenses.