Inverse analysis of laboratory data and observations for evaluation of backward erosion piping process

Sige Peng, John D. Rice

a School of Civil Engineering and Transportation, South China University of Technology, Guangzhou, 510640, China

b Department of Civil and Environmental Engineering, Utah State University, Logan, 84322, UT, USA

Keywords:Backward erosion piping (BEP)Laboratory modeling Inverse analysis Finite element method (FEM)Soil loosening Critical gradient

ABSTRACT An inverse analysis procedure has been developed to interpret collected pore pressure data and observations during backward erosion piping (BEP) initiation and progression in sandy soils. The procedure has been applied to laboratory models designed to mimic the initiation and progression of BEP through a constricted vertical outlet. The inverse analysis uses three-dimensional (3D) finite element method(FEM) to successively produce models of the hydraulic head regime surrounding progressive stages of BEP based on observations at the sample surface and pore pressure measurements obtained from the laboratory models. The inverse analysis results in a series of 3D contour plots that represent the hydraulic-head regime at each stage of the BEP development,allowing for assessing the development of BEP mechanism as well as calculating the critical hydraulic conditions required for various BEP stages to initiate and progress.Interpretation of the results identified four significant stages of the piping process:(1) loosened zone initiation, (2) channel initiation and progression, (3) riser sand fluidization, and (4)loosened zone progression. Interpretation of the hydraulic head contour plots allows assessment of the critical hydraulic gradients needed to initiate and progress various components of the BEP development.

1. Introduction

Backward erosion piping (BEP) is one of the least understood mechanisms of internal erosion,even though it is responsible for a large percentage of the internal erosion failures worldwide(Foster et al.,2000;Vrijling et al.,2010). BEP occurs when non-plastic soil particles are removed at a seepage exit point. The erosion propagates backward toward the seepage source, forming a channel or“pipe” through which soil particles are eroded.

A schematic interpretation of one situation where BEP commonly occurs (through a defect in an overlying lowpermeability layer) is presented in Fig. 1. BEP occurs in this case when flow is concentrated on a defect with an overlying lowpermeability layer. The converging flow creates high hydraulic gradients and flows that enable soil particles to be detached from the underlying sand layer and transported upward through the pipe. BEP may initiate with the formation of a loosened zone that initially forms at the exit location and enlarges with the developing erosion (Fleshman and Rice, 2013, 2014; Allan, 2017). As the hydraulic gradient and flow velocity increase, flow in the defect begins to carry soil particles in suspension and the soil is “fluidized”(Allan, 2017; Van Beek et al., 2013, 2015). With further development,individual soil particles are detached and transported to the exit and the BEP progresses with the formation of one or several pipes. Because the pipe is a path of least resistance, seepage flow from the surrounding soil tends to converge into the pipe. The converging seepage increases both the gradient at the pipe head and the flow within the pipe, enhancing the erosion and transportation potential. With the increasing flow, the pipe may propagate toward the flow source,forming an open pipe and eventually leading to the collapse of the structure.

Assessment of piping potential can date back to the early 1900s.Bligh(1910,1913)and Lane(1935)developed empirical formulae to determine the critical differential head across a structure and presented the empirical resistance for different types of embankment soils.In the early 1900s,Terzaghi and his co-worker(Terzaghi,1922; Terzaghi and Peck, 1948) compared the vertical hydraulic gradient at the ground surface (the exit gradient) with the critical gradient icrneeded to initiate erosion in the affected soil. Terzaghi defined the critical gradient as the ratio of unit weight of soil buoyant (γb) to the unit weight of water (γw), i.e.

Fig. 1. Schematic illustration of BEP and the mechanisms of BEP development with expanded detail of the area represented by the tests in this study.

Although Terzaghi emphasized that his derivation was intended to model the heave process (Terzaghi and Peck,1948), the critical gradient in Eq. (1) has often been used to assess BEP potential.

Sellmeijer and his co-workers from Delft Geotechnics Laboratories and Delft Hydraulics in the Netherlands (e.g. De Wit et al.,1981; Sellmeijer, 1988; Koenders and Sellmeijer, 1992; Weijers and Sellmeijer, 1993; Technical Advisory Committee on Flood Defenses, 1999) performed a wide variety of BEP tests in various flumes to determine the critical seepage condition for BEP development beneath a structure. Several experiments indicated that BEP can initiate at one differential head and then reach an equilibrium state without further BEP development (progression controlled).More increase in differential head is needed for further erosion progress. Equilibrium state could occur at successive increasing differential heads until a critical state was achieved and BEP progressed beneath the entire structure.Schmertmann(2000)also studied BEP progression in a variety of experimental flumes with the aim of analyzing critical heads needed for BEP initiation and progression. One of Schmertmann’s key findings is that the erosion potential is highly dependent on the uniformity coefficient of the sand.

Van Beek et al. (2011, 2015) performed additional flume experiments and modified Sellmeijer’s model using data obtained in the experiments and from previous literature. Van Beek et al. (2011)performed small-, medium- and full-scale experiments with various exit conditions including slope-type and trench-type exit across the entire model (essentially two-dimensional (2D) flow),and a point exit(three-dimensional(3D)flow convergence).It was observed that when BEP is initiated, the differential head with 3D outlet is lower than that with the 2D outlet types, after which the pipe reaches an equilibrium condition and pipe progression is halted(i.e.progression controlled).This behavior is a consequence of flow convergence due to 3D flow. Typically, there is no equilibrium stage in flume experiments with slope or trench-type exits(i.e. initiation controlled).

Allan (2017) performed over 100 flume tests to further investigate the BEP process. One key finding of Allan’s work was the effects of a loosened zone forming prior to the formation of piping channels. Allan used an exponential function to model the increased permeability of the eroding soil near the exit.

Parekh et al.(2016)studied the temporal behavior of backward erosion using the results from large-scale embankment tests at the IJkdijk test facility in the Netherlands,where a differential head was imposed on a full-scale test levee to induce BEP failure. A finite element method (FEM) technique was developed to evaluate pore pressure data from a piezometer array sampled continuously during controlled hydraulic loading. BEP development beneath the embankment was traced by comparing observed pore pressure patterns before and after BEP initiation.

While the researches above have contributed greatly to the understanding of the BEP mechanism in a global scale, further research on mechanism details on the scale of the developing pipe head is still needed.This is especially true where the flow is three-dimensionally constricted into outlets or the head of a developing pipe.The research presented herein was designed to investigate the mechanisms of BEP initiation and early progression under a condition often encountered in the field: exiting into a constricted outlet.In this condition,it is not uncommon to encounter flow that has a large vertical component due to flow converging from depth as well as laterally. This condition is shown in the expanded rectangular area of Fig.1. The vertical component of flow can also be enhanced due to increasing permeability with depth in the sandy layer.

Table 1Properties and characteristics of soils tested.

Fig. 2. Schematic illustration of the testing apparatus.

Fig.3. Locations of pore pressure ports within the sample holder:(a)Top view,and(b)Side view.

This study built upon previous laboratory tests and modeled them using a similar apparatus developed by Fleshman and Rice(2013, 2014) to observe critical gradients with vertical flow through various sandy soils.Using similar methods,Keizer(Keizer,2015;Keizer et al.,2016a,b)tested sandy soils under the condition of inclined exit face and developed relationships to predict critical gradients for BEP initiation under a range of exit face inclinations(Keizer et al., 2016b).

The test in this study investigates the seepage conditions under which BEP will be induced and progress where flow exits into a constricted vertical outlet, mimicking the defect in Fig. 1. An inverse analysis procedure based on FEM modeling is used in this study to analyze the seepage regime surrounding the developing BEP. Input data for analyses include an array of measured pore pressure data within the samples and observations of BEP channel formation at the top of the samples. The resulting 3D models are then used to assess the mechanisms of the initiation and early development of BEP and the critical conditions for BEP further development.

2. Laboratory testing procedure

The laboratory testing procedure was designed to track the initiation and progression of BEP with a constricted entrance through visual observation and an array of closely-placed pore pressure sensors. The data and observations are then interpreted with inverse analysis.

2.1. Materials for testing

Tests were performed on three different sandy soils:(1)graded Ottawa sand(well-rounded silica sand)conforming to ASTMC778-03(ASTM,2017),(2)graded angular sand,angular silica sand with a gradation matching that of the graded Ottawa sand, and (3)uniform, fine-grained (No. 100 sieve) garnet sand. A summary of key properties of the tested soils is presented in Table 1.The specific gravity of garnet sand is much higher than that of the quartz sand(4.05 for garnet versus 2.65 for quartz).

Fig. 4. Test data for graded angular sand: (a) Data of pore pressure and flow plotted versus time; and (b) Normalized pore pressure data plotted versus time.

2.2. Testing apparatus

The laboratory models in this study were conducted using an apparatus similar to the one used in previous studies at the Utah State University (Fleshman, 2012; Fleshman and Rice, 2013, 2014;Keizer, 2015; Keizer et al., 2016a, b) with modification of adding a Plexiglas plate with a round orifice and a riser tube on top of the sample to model a constricted exit.Additional instrumentation was added to capture the complexity of the converging seepage flow.The apparatus imposes an upward gradient on the soil sample to initiate internal erosion, as illustrated in Fig. 2. The soil sample is retained in a Plexiglas sample holder that is sealed in a vertical position between two pressure cells. The system forces upward flow through the soil sample that then converges on the circular orifice at the top of the sample to model the initiation of BEP through a defect in an impermeable clay layer. The head of the high-head reservoir (connected to the bottom pressure cell) is controlled by raising the Mariotte tube while the head of low-head reservoir (connected to the top pressure cell) maintains constant.This system allows for back pressure applied to the system to assist in obtaining sample saturation.

Fig.5. Plot of differential head between sensor PPC verses total differential head across the sample with linear interpolation line and equation for linear (ambient) portion of the plot.

The soil sample is a Plexiglas cylinder (12.7 cm in height and 10.2 cm in diameter) with a plate sealed over the top to form a constricted exit (see Fig. 2). The plate contains a circular orifice at the center with a 5.1 cm high riser tube of the same inside diameter.Tests presented herein were performed using orifice diameters of 1.9 cm (additional orifice diameters are being used for further research). The inside walls of the sample holder were coated with silicone gel for a dual purpose: (1) providing a frictional interface between the soils and the sample holder, and (2) reducing the potential of a preferred seepage path along the sample edges that would occur as a result of larger interstitial voids due to the lack of interlocking with the rigid Plexiglas surface. Two silicone sheets were fixed between the soil and the top plate to improve the contact between the soil and the lid and resemble the sand more closely to clay contact than direct contact with the Plexiglas. The base had screens for retaining the soil particles while allowing water to flow freely into the sample.

Seven pore pressure measurement ports were placed within the sample holder (see Fig. 3). Four ports were embedded within the soil: three located along the vertical axis of the sample at the distance of 1.9 cm(PPA),5.7 cm(PPB),and 9.5 cm(PPC)below the top,and the fourth (PPD)at the distance of 0.95 cm away from the top and shifted 1.14 cm away from the center. Three more ports were located at the top of the soil sample (bottom of the top plate) at radius of 2.86 cm away from center of the riser, radially spaced at interval of 120°.

Differential pressure transducers(Validyne DP15-26)were used to measure differential hydraulic head between the measurement ports and the low-head reservoir (top of sample). A Kobold magnetic flux flowmeter was configured to measure the flowrate between the reservoirs in the tests.The data of pressure transducers and flowmeter were collected by a data logger(Campbell Scientific CR3000, Logan, Utah, USA) and can be observed in real-time on a dedicated computer screen throughout the tests. Video was taken throughout each test and linked to the data using a digital counter in the field of view, thus allowing data correlation with the observed soil behaviors.

2.3. Testing procedure

The tests were conducted with the following procedures:

(1) Soils were placed in the sample holder by dry raining and densified by tapping a metal rod with lift of 1.2 cm.

Fig. 6. Sketches of BEP development from video at key stages of test on graded angular sand (stage letters correlate with Fig. 4).

Fig. 7. A 3D FEM of piping erosion at stage h for graded angular sand.

Table 2Permeability (cm/s) of conducted simulation parameters on graded angular sand.

(2) Two silicone sheets and the lid/riser were fixed on top of the soil sample. The cylinder was over-filled slightly and then densified by compressing the silicon and lid, which was found to be the most effective way to produce a uniform sand density within the soil sample.

(3) The sample was saturated under a partial vacuum by (a)flushing with carbon dioxide, and (b) slowly filling the pressure cells and sample with deaired water from the bottom to the top.

(4) The lead lines of the pore pressure transducer were connected. The saturation process was completed by applying backpressure of 103 kPa(15 psi) to the entire system.

(5) The video recorder and data collector were initiated.

(6) The differential head for each test was set to zero across the soil sample in the beginning and gradually increased by around 1.2 cm until the first movement of sand grains was observed.Then,the increment of loading rate was decreased to about 0.6 cm, allowing the erosion to reach equilibrium under the imposed head at each stage before further increase. The test was continued until the sample completely failed or the sample reached equilibrium with the maximum differential head.

3. Test results and data interpretation

Each test provided a time history of the pore pressure data for each pore pressure port, the total differential head, and the flow data. These data for a test on graded angular sand are plotted in Fig. 4a. In addition, behavior of the developing erosion channels(pipes) observed through the top of the sampler was recorded.Various stages of erosion development (described in Section 3.1)are associated with the data and depicted in Fig. 4a, separated by vertical lines a-n.

Fig. 8. Plot comparing measured pore pressures versus results of FEM models for all sensors at stage h for graded angular sand. (Note: differential head for this stage was 14 cm).

Fig. 9. Average and maximum variations between laboratory testing and FEM results at each stage for three types of sands.

Data interpretation was performed in order to assess the hydraulic conditions (i.e. critical gradients) needed for different stages of BEP development observed during the experiments (i.e.loosened zone initiation,channel initiation and progression,riser sand fluidization,and loosened zone progression).Test data were analyzed by two steps: (1) normalizing the pore pressure data to ambient pore pressure conditions to identify the time when the changes of the flow regime occurred during the tests, and (2)developing 3D models of the seepage regime(i.e.contours of total hydraulic head) using the inverse analysis that utilizes FEM analyses and regression analyses of the results. The total head contours can then be used to assess the critical conditions at the interface between the base soil and the loosened zone, and the channel of various stages of BEP initiation and progression. The data interpretation procedure is presented in Sections 3.1-3.4 by describing the process that is applied to one of the tested soils(graded angular sand).

3.1. Data normalization

Fig.10. Total head contours on the top of the sampler for a test of graded angular sand at various stages(see Figs.4 and 5 for reference of stages):(a)Stage a,ΔH=5.1 cm;(b)Stage c, ΔH = 7.6 cm; (c) Stage d, ΔH = 8.9 cm; and (d) Stage k, ΔH = 17.8 cm.

The term “ambient condition” has been used by Sellmeijer(1988) and Schmertmann (2000) to describe the pore pressure regime that exists prior to any changes in the soil structure caused by internal erosion(i.e.a loosened zone or a pipe).If soil structure does not change,soil permeability within the flow regime will stay constant and the steady-state head variation at any point within the flow regime will vary linearly with the differential heads imposed on the system. This linear behavior allows for easy prediction of the expected hydraulic head values under ambient condition by extrapolating the linear behavior observed at the beginning of the test(before the loosened zone formation changes the flow regime).With the normalized data,it is easier to detect the effects of soil loosening and erosion on the observed behavior.

The raw data collected for the test on graded angular soil over elapsed time are presented on Fig.4a.The data include(1)the total differential head across the sample, (2) the differential head between each pore pressure port and the low-head reservoir,and(3)the measured flowrate.

To quantify the ambient behavior and extrapolate the ambient conditions beyond the end of ambient behavior, equations were developed to describe the linear behavior of the ratio of (1) the differential head between the internal sensors and the top of the sample (Δh1, Δh2,Δh3,ΔhD,ΔhA, ΔhB, and ΔhC) and (2) the overall differential head imposed on the sample,ΔH.For example,the data for sensor PPC versus ΔH are plotted in Fig.5.By regression analysis of the linear (ambient) portion of this plot at small ΔH values(below 4.9 cm),a linear equation describing the expected behavior when no internal erosion occurred (i.e. ambient behavior) was deduced. The linear extrapolation of ambient conditions is presented as the dotted line in Fig. 5 along with the measured data(solid line). Only the portion of the data up to a total differential head of 4.9 cm fitted to the linear equation.The deviation from the linear after that was due to the effects of erosion(i.e.soil loosening and pipe development). Eq. (2) is derived by dividing the differential head between sensor PPC and the low-head reservoir (ΔhC)by the equation for the extrapolated ambient behavior (ΔhE-C).

Fig. 11. Total head contours on a vertical slice through the sampler for a test of graded angular sand at various stages (see Figs. 4 and 5 for reference of stages): (a) Stage a,ΔH = 5.1 cm; (b) Stage c, ΔH = 7.6 cm; (c) Stage d, ΔH = 8.9 cm; and (d) Stage k, ΔH = 17.8 cm.

Eq. (2) is used to calculate the normalized values of the data from PPC and the resulting equation for the normalized differential head, ΔhN-C, for one specific test on graded angular sand. The normalized differential head for PPC, ΔhN-C, versus time is then plotted in Fig. 4b. Similar normalization was performed on the remaining sensors. The resulting values of all seven sensors are plotted in Fig. 4b.

Visual observations of erosion progress were made through the video recorded at the top of the sample top plate. Vertical lines representing key stages (a-n) of BEP progression are plotted in Fig.4a and b and the sketches of the BEP progression for each stage are provided in Fig.6.By comparing the vertical lines in Fig.4 with the sketches in Fig. 6, the effect of BEP progression on the normalized differential head plots in Fig.4b is clearly displayed.For example,between stages b and c,significant pipe growth occurs in the direction of PP2 in Fig. 6. This pipe growth corresponds to a deep drop of ΔhN-2in Fig. 4b, thus illustrating how the pipe formation has caused the reading at PP2 to deviate from the ambient conditions. Similar deviations can be observed between stages e and f and stages g and h due to the growth of pipes in the directions of PP2 and PP1, respectively.

3.2. Inverse analysis methodology

As discussed previously, the BEP process consists of a complex interrelationship among soil loosening, particle detachment, soil transport, and pipe formation. While pipe formation can be observed on the surface of the sand layer,downward loosened zone progression must be inferred by changes in the pore pressure sensor array. That is, further deviations from the ambient conditions which are not induced by the observed pipe formation are due to progression of the loosened zone.To interpret the effects of BEP progression that occurred below the surface, an inverse analysis technique,3D FEM analysis,consisting of iterative assumptions on the extent of soil loosening,was developed to assess the growth of loosened zone during the stage,and back-check the validity of the assumed extent of BEP progression.

The 3D FEM models were developed to model the effects of piping progression at various stages of BEP progression using the computer program SVFlux (Soilvision Systems Ltd., 2009). The boundary conditions for the general FEM model consist of the following aspects: (1) the outlet on the top of the sample was modeled as a zero-pressure boundary;(2)the remainder of the top of the sample was modeled as a no-flow boundary;(3)the bottom was modeled as a constant-head boundary to match the differential head imposed across the sample for each stage;and(4)the sides of the model were no-flow boundaries. The model for stage h with total differential head of 27.9 cm is presented in Fig.7.Considering the developing erosion in successive stages, the models had to be modified. Description of the evolution of model development was presented below. Parameters for various model components of all stages are tabulated in Table 2.

(1) Step 1

A base model, representing the ambient soil conditions, was developed to represent the sample prior to any soil loosening or erosion.The initial void ratio and hydraulic conductivity of the soil at the beginning of the test were measured and confirmed at the early stage of the tests with the flowmeter measurements.The base model consisted of 15,205 nodes and 10,522 elements.However,as changes were made to the model considering the developing erosion,the meshes changed along with these modifications.While soil density and hydraulic conductivity were expected to be uniform throughout the base FEM model,it was necessary to add a thin soil layer with higher hydraulic conductivity on top of the model to match the seepage regime measured during the early (ambient)test stages. This layer was referred to as the adjustment layer and was modeled with a layer with thickness of 0.03 cm along the entire top of the model.The hydraulic conductivity of this layer was adjusted until the model results matched those measured in the early stages of the experiments. The reason that this zone has higher permeability in the physical models is due to a combination of two factors:(1)the contact between soil grains and the rigid top of the sample holder resulting in larger interstitial voids,and(2)a zone looser than the rest of the sample which is induced by unavoidable sample preparation issues (i.e. no confinement in the upper layer as the sample is vibrated). It should be noted that the thickness of 0.03 cm used in the model is arbitrary and it is the transmissivity of the layer (hydraulic conductivity multiplied by height) that is actually being modeled.

(2) Step 2

Fig.12. Schematic illustrations of a scenario where mechanisms of erosion progression are observed in FEM models: (a) loosened zone initiation, (2) channel initiation and progression, (3) riser sand fluidization, and (4) loosened zone progression.

As the differential head increases, changes in soil structure are detected by some of the pore pressure sensors deviating from ambient conditions(see Fig.4b,stage a).Initially,no pipe channels were observed and deviation from ambient conditions was attributed to soil loosening near the exit. Hydraulic conductivity of the loosened-zone was hypothesized to be five times that of the base soil based on previous assessments of loosened zone permeability observed in tests of simpler models without constricted exits (e.g.Fleshman and Rice, 2013, 2014). The lateral extension of the loosened zone was assessed by modeling a cylindrical loosened zone centered below the orifice and adjusting the depth and diameter until the model matched the pore pressure regime measured by the sensors.

(3) Step 3

At later stages,piping channels are observed around the orifice as indicated in Fig. 6, requiring further model modification. The shape of the horizontal channel is based on observations from the video (see Fig. 6). The channels were modeled with a depth of 0.09 cm (about 2.5 times the size of average grain), and the hydraulic conductivity of the channels was adjusted until they corresponded with the experimental results.As with the adjustment layer,the height of 0.09 cm was arbitrary and the transmissivity of the channels was modeled by adjusting the hydraulic conductivity.These back analyses relied heavily on changes occurring in the uppermost three pore pressure ports(PP1,PP2,and PP3).It should also be noted that the transmissivity used in the model represents the total effect of the open pipe channels and loosened zone that may form around the channels.

The hydraulic conductivity values used to model piping channels at various stages in Fig. 6 are presented in Table 2. The transmissivity of the channels changes between the different stages (see Table 2) as the depth of the channels changes with the erosion progress. Other variations in modeled hydraulic conductivity are likely attributable to: (1) clogging of portions of a channel due to deposition or soil loosening from below, (2)variations in soil erodibility along the channel, and (3) local variation in the soils. The largest variation in channel hydraulic conductivity occurred in channel 2 during stages e-g, where it increased to be about 4 times the average value. It was hypothesized that the temporary increase occurred due to the channel expansion in multiple directions during these stages. With expansion occurring in multiple directions, the flow in the pipe center increased, increasing erosion potential and deepening the channel. These resulted in an increasing channel conductance(represented as increasing effective hydraulic conductivity in the FEM model).

Fig.13. Vertical gradients at base of loosened zone during a test for (a) garnet sand; (b) graded Ottawa sand; and (c) graded angular sand.

(4) Step 4

Once the piping channels was modeled based on the upper three pore pressure ports(PP1,PP2,and PP3),any residual changes in the pore pressure regime were attributed to the increases of the loosened zone dimension which was then assessed through further inverse analysis,whereby the dimension of the loosened zone was adjusted to satisfy the measured heads in the lower pressure ports(PPA, PPB, PPC, PPD). In many cases, the results of our analyses indicated that the loosened zone progression was halted once piping channels initiated, and growth restarted only when the channel progressed near the sample edges, a phenomenon discussed later in Section 4.

(5) Step 5

With further piping progression,sand particles began to gather in the riser. This makes it necessary to add a soil zone within the riser to the model effects of the loose and fluidized sand in the riser(see Table 2).The shape and height of the riser sand were based on the pictures and video taken during the test. Originally, hydraulic conductivity of the riser sand represents the effective hydraulic conductivity of fluidized sand boiling in the riser. As the fluidized soil in the riser is not truly Darcian, it is modeled with a high hydraulic conductivity to simulate the added resistance of the fluidized soil.As sand particles started to exit from the riser,the effective hydraulic conductivity increased due to the reduced density of the fluidized mass (fewer particles in suspension).

(6) Step 6

At later stages of the test, the channels stopped growing laterally, yet the volume of riser sand increased and deeper sensors continued to deviate from ambient conditions. Thus, the size of the loosened zone was modified to adjust to these continued changes.

Fig.14. Vertical gradients at base of channels during a test for (a) garnet sand; (b) graded Ottawa sand; and (c) graded angular sand.

3.3. Validation of inverse analysis results

Comparisons were made to assess how well the FEM models developed in the inverse analysis process correlated with actual measured values. Comparison between experimental results and measured ones is presented in Fig. 8 for a test stage on graded angular sand (stage h with a total differential head of 14 cm, see Figs. 4a and 6). Fig. 8 shows that a good agreement between the experimental results recorded by all sensors and FEM results.Similar plots were prepared for all test stages.

A summary of results from plots, similar to Fig. 8, for tests on all three types of soils is presented in Fig. 9, where the average and maximum differences between the experimental and FEM results for all sensors versus the total differential head are plotted. The average deviation was calculated by averaging the percent differences between experimental and FEM results for all sensors at each stage of the test.The maximum deviation was the maximum percent difference for all sensors at each stage. For all stages and tests, the average deviation is less than 4% and the maximum deviation of the seven sensors is no more than 14%.Greater variations occurred at the early test stages and were likely due to the error magnitudes which represented a higher percentage of the total head. Furthermore, the maximum deviations generally came from one of the three sensors at the top of the sample;while the FEM model assumed these three sensors to be identical since piping channels had not been observed yet(this is not necessarily true since differential loosening may have occurred before being observable). Thus, a portion of these deviations can be attributed to random variations within the sample.

Fig.15. Horizontal gradients into sides of channels during a test for (a) garnet sand; (b) graded Ottawa sand; and (c) graded angular sand.

3.4. Use of analysis results

Effects on the pore pressure regime caused by progressing BEP stages and mechanisms (i.e. loosened zone initiation and progression, piping channel initiation and progression, and riser sand initiation and progression)can be observed in the FEM results from the completed inverse analysis procedure.FEM analysis results for four stages of a test on graded angular sand are presented in Figs.10 and 11. Total head contours located right below the sample top plate are presented in Fig. 10, while plots on a slice through the sample center are presented in Fig.11. The first signs of BEP initiation are observed in Figs.10a and 11a representing stage a (total differential head of 5.1 cm), and the extension of a very small loosened zone is identified through inverse analysis. Figs.10b and 11b present the results after two piping channels developed at stage c(total differential head of 7.6 cm).As sand from the channels and loosened zone eroded into the exit, soil particles started entering in the riser at stage d(total differential head of 8.9 cm).The extension of the channels is described in Fig.10c,and changes in the loosened zone and sand appearing in the riser are described in Fig.11c. At stage k (see Figs.10d and 11d, total differential head of 17.8 cm),as fluidized sand started exiting from the top of the riser,the loosened zone again enlarged at a significant rate. The reason accounted for reactivation of loosened zone expansion is believed to be the slowing expansion of pipe channel for reaching the edges of the testing apparatus.

4. Discussion of results

Visual observations and results of inverse analysis described above indicate that the BEP progression observed in the laboratory tests consists of four observable stages: (1) loosened zone initiation, (2) channel initiation and progression, (3) riser sand fluidization, and (4) loosened zone progression. Schematic of interaction among these stages of BEP progression in the laboratory models is provided in Fig.12.

In Fig.12a,a small loosened soil zone is formed below the exit due to flow concentrated on this area. The inverse analysis of the pore pressure data indicated that the loosened zone gradually enlarged through the early stages of hydraulic loading until horizontal gradients in the upper sample became large enough to initiate the channel formation. Initial development of piping channels is illustrated in Fig. 12b. As channels formed, seepage flow that had been concentrating in the loosened zone converged into the high-conductivity channels, thus decreasing hydraulic gradients at the loosened zone boundary and slowing or halting its enlargement. With increased total differential head, the channels propagated with little or no enlargement of loosened zone and soil particles began to enter the riser due to the soil particles eroding from the channels (see Fig.12c). The riser sand was assumed to have a constant hydraulic conductivity at the beginning as the sand rises when its density decreased due to the removal of particles. It should be noted that the hydraulic conductivity used for the fluidized soil in the riser is not equal to that for Darcy flow since it is modeling the effects of a dense fluid on top of the soil rather than a true Darcian seepage. Nonetheless,assigning a hydraulic conductivity value to this layer has the same effect in the FEM analyses as the dense fluid does. Once the channel development reaches the edge of the sample holder, the loosened zone again begins to increase in width and depth again,as illustrated in Fig.12d.

The inverse analysis results presented herein allowed observation of the hydraulic regimes during development of BEP,including hydraulic gradients occurring at the interfaces between the various components.Based on such observations,assessment of the critical hydraulic condition, under which the erosion was initiated and progressed, could be conducted. The vertical gradients at the loosened zone for tests of three types of soils are plotted in Fig.13,along with annotation of piping progression.These plots exhibited similarities and differences of their behaviors. The general increasing trend of gradient at the loosened zone was interpreted to be induced by the increasing effective stress at this zone base as a result of the increasing depth and effective stress.In the latter two plots(see Fig.13b,c),the gradient at the loosened zone base actually decreases for a portion of the test. The reason for this is that seepage flow converges into the expanding piping channels instead of being concentrated in the loosened zone.

The channel initiation and progression mechanisms are also important for understanding the development of BEP. Vertical gradients at the center of each channel in the three tested soils are plotted in Fig. 14. Unlike gradients into the loosened zone base,gradients into the channel base stayed relatively constant throughout the tests.This is likely due to an equilibrium state along the channel base, allowing the channel size to increase with available seepage flow, thus keeping the gradient relatively constant.

Plots of the maximum horizontal gradients at the edge of developing pipe channels for the three types of soils are presented in Fig.15.Gradients were calculated based on hydraulic head drops occurring at the distances of 0.3 cm and 0.5 cm away from the pipe channel edge.As shown in Fig.15,a moderate amount of horizontal gradient variation exists for each soil type. This variation is likely due to a number of factors including: (1) local soil structure and density variation, (2) pipe channel configuration and associated concentration of flow, and (3) interaction with adjacent BEP components development. Assessment of how these factors affect the critical gradient for pipe development would require additional tests and analyses.

Nevertheless,some complex behaviors such as critical gradients at pipe channel heads and behaviors in risers need further studies,which are beyond the scope of this paper.

5. Conclusions

This paper described an inverse analysis technique used to interpret observed behaviors and pore pressure data collected during laboratory tests. We studied the mechanisms of BEP initiation and early progression in sandy soils when a defect in an overlying soil blanket was reported. The flow entered the model vertically at the base and then converged into the soil before exiting through the constricted exit.This is intended to roughly model the converging conditions,although some amount of horizontal flow is likely in most field conditions when it enters the model area. The purpose of the study was to demonstrate an inverse analysis technique for processing the test data,and to provide insights into the complex mechanisms of BEP.

The observations and pore pressure measurements were interpreted by rectifying the observed progression and changes in the flow regime using 3D FEM analyses.The inverse analysis technique was shown to produce 3D models that accurately matched the observed and measured conditions at various stages of BEP development. From such models, critical hydraulic conditions and hydraulic parameters needed to initiate and progress at various stages of BEP were obtained.

Four significant BEP development stages were identified through the comparison between laboratory tests and FEM model:(1) loosened zone initiation, (2) channel initiation and development,(3)riser sand fluidization,and(4)loosened zone progression.The four stages represent several mechanisms that occur as the soil and flow regimes react to the increasing gradients. Stages 1 and 4 represent the initiation and extension of loosened zone in the samples. Stages 2 and 3 illustrate how flow concentration is defused from the loosened zone to the progressing channels,resulting in slowing or halting of the loosened zone expansion.

The observations of the BEP initiation and progression process provide deeper understanding of the BEP process. Further researches with a wider variety of soil types, orifice diameters and riser heights are needed to provide more insights into the mechanisms of BEP and the capability to predict the occurrence of BEP.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The first author would like to acknowledge the support from the South China University of Technology for the PhD short-term visiting project.

List of symbols

CuCoefficient of uniformity

D50Cumulative particle size grain size

GSSpecific gravity

icrCritical gradient needed to initiate erosion in the affected soil

Δh1Differential head between sensor PP1 and the low-head reservoir

Δh2Differential head between sensor PP2 and the low-head reservoir

Δh3Differential head between sensor PP3 and the low-head reservoir

ΔhDDifferential head between sensor PPD and the low-head reservoir

ΔhADifferential head between sensor PPA and the low-head reservoir

ΔhBDifferential head between sensor PPB and the low-head reservoir

ΔhCDifferential head between sensor PPC and the low-head reservoir

ΔH Total differential head across the sample

ΔhN-CNormalized differential head

φiInternal angle of repose

Ψ Sphericity

θ Roundness