3D finite element method (FEM) simulation of groundwater flow during backward erosion piping

Abstract

Backward erosion piping is an important failure mechanism for cohesive water retaining structures which are founded on a sandy aquifer. At present, the prediction models for safety assessment are often based on 2D assumptions. In this work, a 3D numerical approach of the groundwater flow leading to the erosion mechanism of backward erosion piping is presented and discussed. Comparison of the 2D and 3D numerical results explicitly demonstrates the inherent 3D nature of the piping phenomenon. In addition, the influence of the seepage length is investigated and discussed for both piping initiation and piping progression. The results clearly indicate the superiority of the presented 3D numerical model compared to the established 2D approach. Moreover, the 3D numerical results enable a better understanding of the complex physical mechanism involved in backward erosion piping and thus can lead to a significant improvement in the safety assessment of water retaining structures.

Keywords

backward erosion piping / groundwater flow / 3D finite element method (FEM)

Introduction

A sandy aquifer situated below a cohesive water retaining structure leads to a concentration of groundwater flow toward a downstream exit opening. This may cause an emission of sand grains at that particular location. The erosion of sand grains results in the formation of shallow pipes in the sand layer right below the water retaining structure. These pipes do not collapse under gravity because of the bridging nature of the cohesive material. The term ‘backward erosion piping’ is designated to the growth direction of the pipes which is opposite to the flow direction, i.e., from downstream to upstream. Eventually, the pipe will form a direct connection between upstream and downstream, which obviously leads to a facilitated water transport and to the action of accelerated erosion. The pipe dimensions finally reach an unbridgeable size for the cohesive water retaining structure, resulting in a (partial) collapse.

Different models exist for the safety assessment of a levee regarding piping failure. Nowadays, Sellmeijer’s prediction model is often used to estimate the critical gradient at which backward erosion piping leads to dike failure, see Eq. (1) [1-3].

(1){\displaystyle \frac{H_{\mathrm{c}\mathrm{r}}}{L}}=F_{\mathrm{R}}{F}_{\mathrm{S}}{F}_{\mathrm{G}},

with F_R, F_S and F_G the resistance factor, the scale factor and the geometrical shape factor respectively:

(2)F_{\mathrm{R}}=\eta {\displaystyle \frac{{\gamma }_{\mathrm{p}}^{\prime }}{{\gamma }_{\mathrm{w}}}}\tan \theta ,

(3){\left({\displaystyle \frac{RD}{R{D}_{\mathrm{m}}}}\right)}^{\mathrm{0.35}}{\left({\displaystyle \frac{U}{U_{\mathrm{m}}}}\right)}^{\mathrm{0.13}}{\left({\displaystyle \frac{KAS}{KA{S}_{\mathrm{m}}}}\right)}^{-\mathrm{0.02}}{F}_{\mathrm{S}}={\displaystyle \frac{d_{\mathrm{70}}}{\sqrt[\mathrm{3}]{\kappa L}}}{\left({\displaystyle \frac{d_{\mathrm{70}\mathrm{m}}}{d_{\mathrm{70}}}}\right)}^{\mathrm{0.6}},

(4)F_{\mathrm{G}}=\mathrm{0.91}{\left({\displaystyle \frac{D}{L}}\right)}^{\frac{\mathrm{0.28}}{{(D/L)}^{\mathrm{0.28}}-\mathrm{1}}+\mathrm{0.04}}

and the submerged unit weight of particles, the unit weight of water, η the coefficient of White, θ the bedding angle, d₇₀ the particle diameter, D the height of the sand layer, L the width of the structure, RD the relative density, U the uniformity, KAS the roundness, and κ the intrinsic permeability. The formula was obtained by a combination of analytical formulae, 2D FEM simulations for modeling groundwater flow, observations from piping experiments on critical pipe lengths and curve fitting. The criterion for onset of erosion is the equilibrium of sand grains at the pipe bottom under the action of water flow within the pipe according to White [4] in 2D. In situations which are similar to those forming part of the study that led to the formula, the prediction model is successful. However, recent experiments indicate that this original Sellmeijer formula fails to predict the critical head correctly under different constraints and an empirical correction has been proposed [3]. A possible theoretical origin for these poor predictions could be the assumption that backward erosion piping can be described as a 2D phenomenon.

In this study, the consequences of the 2D approach for the piping phenomenon are investigated by means of a detailed comparison between 2D and 3D finite element method (FEM) simulations of the groundwater flow. In addition, the influence of the model width on the exit velocity is investigated in order to determine the area of influence actively participating to the piping process.

2D approach for piping initiation

Basically, the 2D assumption invokes that the water, which actively contributes to the erosion mechanism, originates from a vertical 2D plane (see Fig. 1(a)).

In our numerical setup (Abaqus 6.12), the sand is defined as a porous material with a specific permeability, while a hydraulic gradient is applied by means of pore pressure boundaries (both upstream and downstream). 4-node plane strain quadrilateral elements (CPE4P) are employed with varying sizes across the model in order to achieve accurate results, while keeping the computation time acceptable. The erosion length E, the dimension of the downstream exit d, the height of the sand layer H and the hydraulic head drop ∆H are respectively 340, 6, 100 and 20 mm. These dimensions have been specifically chosen in order to compare with the performed small-scale experiments. The permeability of the used sand type in these experiments (denoted as Eastern River sand) is 2.4 × 10^-4 m/s.

A contour plot of the flow velocity is shown in Fig. 2. The water flow toward the open exit can be readily observed. Piping initiation is triggered once the vertical exit velocity (VE) in the outlet opening (shown in Fig. 3) exceeds a particular threshold over a certain distance [5]. Note that for infinitely small mesh elements, the maximum velocity would approach infinity at the boundaries, as prescribed by theory. Here, the mesh elements are chosen small enough, such that the exit velocity at a distance of 0.3 mm from the hole boundaries does not change for finer elements.

3D approach for piping initiation

3D Numerical simulation of a ditch type exit

As a first step toward 3D simulations, the ditch-type exit (see Fig. 1(b)) has been simulated. The CPE4P elements have been replaced by their 3D counterparts, namely the C3D8P elements. In principle, this 3D simulation should yield the same results as obtained for the 2D approach because of translational symmetry, at least when assuming an infinite width of the ditch. The latter has been ensured by application of suitable symmetry conditions. The result for the computed vertical exit velocity is added to Fig. 3. It can be readily verified that similar results as in the 2D simulation are obtained. Hence, these results validate the implementation of the 3D FEM.

3D Numerical simulation of hole type exit

In reality, the exit type for backward erosion piping is typically a hole. Therefore, an impermeable downstream layer with a hole type exit with dia-meter d = 6 mm (see Fig. 4) is considered here. Similar as for the 2D simulation, the erosion length E, height H and head difference ∆H are respectively 340, 100 and 20 mm. The width of the 3D model is set to 30 cm, which equals the dimension of the performed small-scale experiments.

The contour plot of the flow velocity in the xz-plane is shown in Fig. 5(a). Although the global flow velocity within the sand is smaller compared to the 2D situation, the overall exit velocity is significantly higher (see Fig. 3 for a comparison of the flow velocity profiles). This can be understood in terms of water supply along the 3rd dimension of the model, as illustrated in Figs. 6(a) and (b).

Figure 5(c) shows contour plots of the flow velocity in top views at different z-positions. These results clearly show that the initial uniform inflow of water upstream is further downstream concentrated toward the hole. Likewise, the contour plot in the yz-plane (Fig. 5(b)) illustrates the confluence of water toward the exit opening. The 3D results indicate that the exit velocity is much higher compared to the 2D results, see also Figs. 2 and 3. Hence it is clear that a safety assessment for piping initiation, based on a 2D approach, will result in inaccurate results for the start of the pipe and consequently in an unsafe levee design for the natural case of an exit hole.

To further investigate the exact influence of the 3rd dimension on the computed exit velocity, several 3D simulations are performed in which the width of the model has been varied. Some of these results are displayed in Fig. 7(a). It can be observed that the exit velocity VE significantly increases for an increasing width W. The rising trend becomes less pronounced after a certain threshold. Hence these results suggest that the horizontal water flow, rather than the vertical water flow, is the driving force to initiate the 3D phenomenon of piping. For brevity, we suffice by noting that equivalent results are obtained for different erosion lengths. Even more, it was found that the value of the particular threshold relates to the erosion length. This is shown in Fig. 7(b), in which an erosion length of 34 mm is considered.

It can be stated that a width W equal to the erosion length is sufficient but also necessary to accurately capture the 3D phenomena associated with piping initiation. Hence the width of a small-scale experimental setup for research purposes on piping initiation should be carefully chosen.

The previous results might explain the periodically spaced outflow openings which were observed in nature [6] in the absence of an impermeable top layer downstream of the levee: when the distance along the dike becomes too large, the water is not sufficiently attracted by an existing outflow opening (see Fig. 6(b)), and thus creates a new outflow opening at a distance from the original hole equal to the influence width. Conversely, multiple outflow openings at distances smaller than the influence width are less susceptible to piping because they have to share the water influx, resulting in smaller flow velocities (see Fig.8).

Experimental validation of the 3D approach for piping initiation

Backward erosion piping has been experimentally investigated [5,7] by means of a small-scale physical model. The experimentally measured flow amounts 1.55 mL/min, while the numerically calculated flow amounts 1.41 mL/min. Note that the simulation has been performed with the same dimensions and similar parameters as the small-scale experiment (∆H = 20 mm). The small deviation between experiment and FEM is probably attributable to small irregularities and uncertainties in the experimental setup.

3D approach for piping progression

The importance of the 3rd dimension on the initiation of backward erosion piping has been discussed in the previous paragraphs. The question arises whether its influence is similar for the further development of the pipe.A variety of pipe development stages (Table 1) which evolved during experimenting are considered. The pipe has been modeled as a simplified straight canal, having a fixed triangular cross-section.Its volume and dimensions are estimated based on the volume of the deposited sand, which actually forms a crater, the pipe length and horizontal area which are measured based on photos of the experiment (see Fig. 9).

As a first approach, the pipe has been modeled as a porous material with a high permeability. For any given pipe permeability however, poor agreement with the experiment is obtained for all 7 pipe development stages (Table 1). A reasonable alternative was obtained assuming an increased sand permeability in the area surrounding the pipe. As described in Ref. [8], sheet flow occurs at the interface between sand and water (bottom of the pipe) and thus provokes a higher porosity of the sand which results in an increased permeability. Moreover, semi 3D experiments in which the piping process is visualized from the side [9], reveal a slantwise rolling and falling of the detached sand grains on the pipe wall and bottom. This, combined with the upward force on the grains caused by the gradient might cause a less dense packing at those locations. For a certain set of parameters (permeability of the pipe= 0.1 m/s, permeability of the sand with increased porosity= 0.01 m/s, sand with increased porosity: depth= 15 mm and width between 17 and 24 mm around the actual pipe), the assumption of a more permeable sand layer around the pipe yields the results in Fig. 10, showing the flow as a function of the pipe volume (representing different development stages) with a full line for the numerical results and a dashed line for the experimental measurements. The agreement between the numerical and the experimental flow is relatively good for all 7 pipe development stages, especially when bearing in mind the simplifications that are made on the geometry of the pipe. It is acknowledged however, that the numerical model could be optimized with more extensive experimental and numerical research to fully capture the mechanisms involved.

The influence of the model width for the differentpipe development stages becomes clear by plotting the exit velocity at the center {VE}_{3D}^{c\mathrm{ent}er} divided by the maximum exit velocity for the standard 2D situation VE_c,2D (see Fig. 11). It can be observed that two trends are present. After initiation of the pipe (stage 1), the exit velocity increases for a width W which extends the erosion length (stage 2 to 5). This can be understood by taking into account that a pipe transports more water which must be supplied from a larger area of influence (see Figs. 6(b) and (c)). When the pipe has reached the upstream side (stage 6 and 7) however, the presence of large inflow areas becomes less significant since the contribution of water supply through the sand is insignificant with respect to the direct water supply through the pipe itself (see Fig. 6(d)).

These results clearly suggest that the 3D nature of the water flow becomes even more significant during the propagation of a pipe. From Fig. 11, it can be concluded that the width W of the setup should be at least 3 times the erosion length in order to accurately examine the growth of a pipe, at least for the assumptions made in this study.

The previous reveals that not only piping initiation, but also piping progression is a full 3D phenomenon.

Discussion

The presented results demonstrate that treating backward erosion piping in small-scale tests as a 2D phenomenon is an error-prone assumption. However, this conclusion is not “a characteristic of small-scale testing”, but it is a principle that should be extended to full-scale situations and real dikes. Moreover, this conclusion suggests that the already very advanced current design formulas that take into account numerous parameters miss an important aspect: the contribution of water flow in the third dimension and the consequent impact on erosion.

It should be noted that the results of this study are not directly applicable to the design calculations of dikes. More research (such as the influence of physical boundary conditions and the criterion for the onset of erosion during the progression phase) needs to be done to fill the existing gap of knowledge on backward erosion piping before the current design formulas can be improved. We believe that the results of this paper demonstrate the importance of treating backward erosion piping as a 3D phenomenon when conducting further research to eventually develop a new formula which is of practical use.

Conclusions

A numerical approach of the groundwater flow leading to the erosion mechanism of backward erosion piping has been presented and discussed. By means of both 2D and 3D FEM simulations, the 3D nature of backward erosion piping has been demonstrated. Large 3D volumes are drained by the pipe, resulting in a significantly larger exit velocity, which in turn holds a greater risk of erosion. Therefore, 2D prediction models may be unreliable.

Furthermore, the 3D computations indicate that a single outflow opening within an impermeable layer is more susceptible to the occurrence of piping than multiple outlet openings or a ditch type exit. In case of a ditch type exit or in the absence of an impermeable top layer downstream of the levee, the numerical results reveal that periodically spaced holes are formed (and thus the 2D situations evolves toward a 3D situation)in order to take full advantage of the available water.

The results of the presented 3D numerical investigation also provide a minimum bound for the width of the small-scale experiments in order to accurately study both piping initiation and piping progression. The former requires a minimum width equal to the erosion length, while the latter requires a minimum width of three times the erosion length.

References

[1] SellmeijerJ B. On the mechanism of piping under impervious structures. Dissertation for the Doctoral Degree. Delft: TU Delft, 1988

[2] SellmeijerJ B, KoendersM A. A mathematical model for piping. Applied Mathematical Modelling, 1991, 15(11–12): 646–651

[3] SellmeijerH, de la CruzJ L, van BeekV M, KnoeffH. Fine-tuning of the backward erosion piping model through small-scale, medium-scale and IJkdijk experiments. European Journal of Environmental and Civil Engineering, 2011, 15(8): 1139–1154

[4] WhiteC M. The equilibrium of grains on the bed of a stream. Proceedings of the Royal Society of London, 1940, 174(958): 322–338

[5] Van BeekV M, van MeertenH, NugrohoD, VandenboerK. SBW Piping 5A: Model Development for Initiation of Piping. Deltares, 2012

[6] MieselD. Rückschreitende Erosion unter bindiger Deckschicht. Berlin, 1978

[7] Van BeekV M, KnoeffH, SellmeijerH. Observations on the process of backward erosion piping in small-, medium- and full-scale experiments. European Journal of Environmental and Civil Engineering, 2011, 15(8): 1115–1137

[8] ShamyEl, AbdelhamidU, KruegerP, AnZ. A Particle-based Model of Flow-Induced Scour. In ICSE6. Paris, 2012, 8

[9] Van BeekV M, van EssenH, VandenboerK. SW Piping 5B: Progression of Piping. Deltares, 2013

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg