SIMULATING LEVEE EROSION WITH PHYSICAL MODELING VALIDATION

This paper studies rill and gully initiation and propagation on levees, dams, and general earth embankments. It specifically studies where these erosion features occur, and how long a particular embankment can sustain overtopping before breaching and catastrophic failure. This contrasts to previous levee erosion analysis, which has primarily concerned the final effects of erosion, such as soil loss, depth of scour and breach width. This paper describes the construction of scaled-down physical models of levees composed of different homogeneous sands, as well as sand- clay mixtures, and their laboratory testing. A 3-D laser range scanner captured the surface features of the physical model, before and after erosion. The resulting data is utilized in developing digital simulations of the rill erosion process. Those simulations combine 3-D Navier-Stokes fluid simulations and a segmented height field data structure to produce an accurate portrayal of the erosive processes, which will be validated by physical modeling.


INTRODUCTION
Levee failures that have occurred as a result of storm surges and flooding events have been primarily due to overtopping, although failure from seepage is also a possible failure mechanism. In either instance, the erosive processes can eventually lead to breaching of the levee and catastrophic damage on the adjacent floodplain. There have been many cases of earth embankment fai lures, for example, Hurricane Katrina in 2005 , where breaching occurred and devastated the surrounding population. Levee failures are preventable, and a better understanding of the ways in which these embankments are designed and fail , so as protect against future failures, is a goal of this research.
The erosion processes described in this paper refer to hydraulic erosion. Small-scale erosion on earthen embankments is being studied, modeled and eventually simulated, with respect to the formation of rills and gullies. Validation of the simulation is a primary focus in this research, so scaled-down model levees are used to perform erosion experiments at I g and later at higher levels of g in a geotechnical centrifuge.
The results of experiments to date, are presented in this paper. Completed testing has been performed at I g using a homogeneous laboratory sand and Nevada sand -kaolin clay mixes. The physical models serve as the basis for developing accurate, digital simulations of the embankment erosion processes. Eventually, different types of soils and soil mixtures will be tested and complex geometries and boundary conditions utilized to quantitatively assess the effects of differing conditions.

RELATED RESEARCH
There is a considerable amount of information pertaining to erosion on earth structures such as levees, dams and embankments both from a civil engineering as well as a computer graphics perspective. Current research on the topic of erosion in the field of civil engineering is primarily associated with developing models that predict final erosive measures (i .e. scour depth, final breach width, total soil loss, etc.) . In the computer graphics field, multiple attempts have been made to simulate hydraulic erosion, chiefly to generate realistic-looking terrain and surface deformation animations as a result of fluid flow . While the research in both fields is beneficial and relevant, neither model the erosion, sediment transport or deposition processes with real physical accuracy capable of predicting the extent of erosion or possible water inundation as a result of breaching.

Erodibility
The erodibility of a soil relates the velocity of the water flowing over the soil to the corresponding erosion rate experienced by the soil. A soil's erodibility was defined as a way to describe the behavior of a soil under erosion conditions. Wan and Fell (2004) describe the development of two erosion rate tests, the Hole Erosion Test (HET) and Soil Erosion Test (SET), which measure a soil's erodibility. Using an Erosion Function Apparatus (EF A), Briaud and his colleagues investigated the erodibility of several different types of soil. The soils were classified into different categories of erodibility based on degree of compaction, erosion rate, water velocity and hydraulic shear stress (Briaud, et al. 2008) . Xu and Zhang (2009) found that in addition to soil type, the degree of compaction plays an important role in erodibility on embankments. The erosion resistance increases with compaction effort, particularly with fine soils. Briaud et al. (2008) deviate slightly from the broad definition of soil erodibility in order to produce a more technically correct definition of the parameter. Since the velocity of the water at the soil-water interface is zero, yet soil is still eroded, soil erodibility is actually based on the hydraulic shear stress. The shear stress changes with water velocity so that it can be defined along the soil-water boundary, incorporating the soil features as well as the water properties along the flow field. This model is ideal for small scale erosion simulation, as it allows for a parameter to be applied to the soil as a field over the entire embankment.

Rill Initiation
Multiple factors influence erosion on an embankment including embankment configuration, flow velocity, slope discontinuities, presence of tailwater, and flow concentration at low points (Pow ledge, et al. 1989). These factors, if present, can impact the formation of rills and gullies on a slope, as well as the shape and speed in which the rill or gully propagates. On a slope, the overland flow first arrives as "sheet" erosion, then causes rill erosion with increasing flux (An and Liu 2009). Bryan and Rockwell (1998) studied agricultural sites near Toronto, Canada and found that significant rill incision typically occurred in early spring, immediately following snowmelt. This relates to the study of levees or earth dams that are adjacent to water bodies. They are saturated or can become nearly saturated rapidly, thereby creating rill initiation conditions.
Rills and gullies will fonn in areas of depression, or in areas where the soil does not have enough cohesion or shear strength to resist the hydraulic stresses from the flowing water. Factors affecting rill characteristics include the stress caused by the flow, roughness of the soil surface, slope gradient and soil erodibility (Mancilla, et al 2005). It was concluded that the most critical determinant of rill development is not threshold hydraulic conditions associated with intense runoff on steep slopes, or areas of depression, but impermeable subsoils that allow surface soils to become saturated (Bryan and Rockwell 1998).

Rill Propagation
After a rill has been initiated in an embankment slope, the initial rill will transport the majority of water and sediment. Occasional tributary rills may form temporarily that supply the main rill with water and sediment, but will taper off as the erosion process continues in the initial main rill (Mancilla, et al 2005). Erodibility within a rill may vary with depth, which can decrease the erosion process in granular soils, as a result of a reduced slope gradient. If a more erodible soil underlies the surface soil, however, the erosion rate in a rill or gully will actually be accelerated (Govers, et al. 2007). Briaud et al. (2008) performed several tests that indicated that the rill erosion occurs first on the land side of the overtopped levee and progressively recedes, leading to eventual breaching. The quantity of soil eroded rapidly increases with the slope gradient, then decreases suggesting a critical slope gradient. If the slope of an embankment has not exceeded the critical gradient, interrill erosion occurs and transports sediment between rills or gullies. The majority of sediment carried by interrill erosion concentrates around rill heads, leading to an increased erosion rate and wider rill width in that area of the rill. Although the incidence of interrill erosion is larger than that of rill erosion, rill erosion is the dominant process on embankment slopes because it is significantly more intense (An and Liu 2009).

Performance Differences Amongst Various Soils
Post Hurricane Katrina field surveys showed that in general, rolled compacted clay fill levees performed well with minor erosion occurring when overtopped, whereas hydraulic filled levees with significant amounts of silt and sand performed poorly. Using good clayey material often required long haul distances that slowed construction progress. So nearby granular material was used instead (Sills, et al. 2008). In cohesive embankments, breaching occurs as a result of headcutting, whereas in granular embankments, surface slips occur rapidly due to seepage on the downstream slope (Xu and Zhang 2009).
Threshold hydraulic stress values tend to be higher on freely drained material. After formation of the water table, however, this value drops, thereby making freely draining granular soil much more erodible (Bryan and Rockwell 1998). Dealing with waste embankments, research by Thornton and Abt (2009) showed that lower clay contents correlated with greater potential susceptibility to the gully erosive process.
Cohesive soils are more resistant to erosion due to high clay content. However, care must be taken when specifying and inspecting the type of clay used. Dispersive clays are an exception because the clay particles spontaneously detach from one another under saturated conditions (Torres 2008). Rockfill and clay embankments are considered to have medium to low erodibility, while silt and sand are considered to have high to medium erodibility according to Briaud' s erodibility classification (Xu and Zhang 2009).

Physically-Based Erosion Simulation
Hydraulic erosion has been accepted as the single most important process in the shaping and development of terrain . Because of this, hydraulic erosion research in computer graphics has focused mainly on terrain generation and animation. The height field erosion simulation performed by Musgrave et al. (1989) and the sedimentation process in the work by Chiba et al. (1998) are examples of terrain generation in computer graphics. In each example, an erosion process is simulated on a terrain to morph and mold it to be more realistic-looking.
Erosion simulations require efficient algorithms that can be run on dynamic data structures in order to capture the small-scale complexity of the process. There are three primary data structures that are often used for erosion simulation: height fields , voxel grids, and layered height fields. Stuetzle, et al. (2010) presented an extension to the layered height field, called a Segmented Height Field (SHF).
Although much work has been done to simulate erosion, very little of it has presented validation of results. Validation of our computer simulation by laboratory experimentation is a primary objective of this research. To our knowledge, validation of computer simulations has not yet been accomplished, though it has been attempted with some success by the Soil Degradation Assessment (SoDA) project (Valette et al. 2006).

PROCEDURES Physical Modeling
Initial overtopping tests were conducted on a half-levee model within an open aluminum box of dimensions 0.356m x 0.61 m x 0.9l4m (14" x 24" x 36"). A piece of plywood, O.013m (W') in thickness, was cut to the dimensions 0.152m x 0.61 m (6" x 24") and sealed in the aluminum box using silicone, which partitioned the space within the box into two distinct zones. The smaller zone measured 0.152m x 0.61 m x 0.216m (6" x 24" x 8 Y2") and was used as a reservoir for the water to rise in and eventually overtop the model levee, which was constructed in the second, larger zone using a moist, medium-well graded laboratory sand having a dry unit weight of 100 pcf and an internal friction angle, <p = 39 .6°.
Constructed in lifts and compacted using a 0.102m x 0.102m (4" x 4") wooden hand tamp, a level base layer 0.076m (3 ") thick was placed first, followed by the half-levee with a 0.127m (5 ") wide crown and 5H: 1 V slope. This slope inclination as per the U.S. Army Corps of Engineers (USACE) Levee Design Manual to prevent damage from seepage exiting the slope on the land side for sand levees (USACE 2000). The elevation of the crown was 0.152m (6"), the same as the elevation of the plywood partition. The final configuration of the model levee in the box left a 0.127m x 0.61m (5" x 24") area, in plan, at the toe of the levee slope that acted as a floodplain. Figure I shows the physical experiment setup. A small aquarium pump was placed 0.0 13m (W') above the floodplain at the farthest point downstream in the box to pump out flood water and allow the overtopping to continue for a longer duration. In both (a) and (b) , the water source is located on the left side of the box, and the pump to remove excess flood water is represented as the black object on the right side of the box.
Overtopping tests were also conducted on a full -levee model in the same aluminum box described above, using the same laboratory sand as well as sand-clay mixes. The sand-clay mixes had a dry unit weight of96 pef. The 0.152m x 0.61m (6" x 24") plywood was replaced with a 0.076m x 0.61 m (3" x 24") piece of plywood that was sealed in the box with silicone, partitioning the box into halves for the laboratory sand testing only. A core was not used during the sand-clay mixture tests, as the levee slopes were flat enough and the soil had sufficient cohesion to prevent seepage failures. The role of the plywood in this setup was to serve as a low-permeability core for the levee. A base layer 0.038m (1 Y2") thick was constructed in lifts and compacted with the hand tamp on each side of the partition. The full-levee was constructed with a 0.203 m (8") wide crown and 5H: I V slopes, so that the elevation of the crown was O.013m (W' ) above the elevation of the plywood core, allowing for some breaching to occur. The pump was placed at the most downstream point, 0.013m ( W') above the floodplain, and water was allowed to overtop the levee for tests with flow rates ranging from 0.0063 Llsec to tests with flow rates of 0.040 Llsec. The majority of tests conducted have used a constant flow rate between 0.0 I 0 Llsec and 0.015 Llsec, however, more extreme cases (very slow flow rate and very rapid flow rate) were investigated to observe the influence of flow rate. In each test, evidence of rill erosion was carefully monitored.

Data Collection
3-D scans were taken before erosion simulation began and immediately following erosion of the levee model. Each scan represented a different terrain elevation of the levee. The scans were taken using a 3-D laser range scanner (Fig. 2a), which provided a point cloud ( Fig. 3c and 3d) with color information for each data point. This data was then processed to ready it for adaptation to the Segmented Height Field (SHF) . To minimize holes in the data from occlusion, two scans were made of each elevation and registered using keypoints on the rigid aluminum box (Smith, et al. 2008). Once each scan for a single elevation was aligned to the same coordinate system, points were then discretized by superimposing a 2-D horizontal grid over the point cloud and snapping each point to its nearest grid space. The data in each grid space was then averaged to create a single height value per grid space, creating a height field. Grid cells not assigned a height value from the scan are interpolated using ODETLAP (Stookey, et al. 2008). This procedure was repeated for each layer, and used to create a single SHF from the layers of height fields, shown in Figure 2.
The renderings in Figures 2b, 3c and 3d were created using the program Mathematica and are colored according to elevation. The relatively equal spacing of the elevation contours in Figure 3c compared to the displaced and irregular elevation contours in Figure 3d illustrate the movement of sediment and the presence of a rill on the slope after water had overtopped the model levee. A slight curve in the elevation contours in Figure 3c is seen on the lands ide slope. The rill formed in that area where the elevation changed slightly, as expected according to previously published findings . Figures 2b and 3d show the rill regression on the levee crown stopping half way across the crown. The receding channel stops at this location because the plywood core had been reached from displacement of sediment and illustrates the effectiveness of a core in at least slowing down the process of a full breach occurring. /.-

Interpolation and Visualization
In order to perform an erosion simulation, a surface must be extracted from the data structure. To do this, the data is converted into a tetrahedral mesh, allowing surface information to be extracted, such as slope, as well as generate surface normals for visualization. These not only improve the quality of the resulting visualization, but also yield more accurate physical simulations by allowing water to flow smoothly down the levee's slopes and through channels cut within the soil.

Half-Levee Setup
Water began to overtop the reservoir and flow through the soil (groundwater, seepage, etc.), thereby saturating the soil, and slowly flowed over the crown of the embankment on the surface. Surface tension was evi dent, as the water on the surface of the crown had boundaries (i.e. the water did not come over the top in one big sheet of water). Once the water had crossed the crown on the surface, rill initiation at the top of the embankment was observed, beginning at the crest of the slope (the intersection of the crown edge and edge of the slope) and eroding its way to the toe of the slope. This formed the primary rill on the slope and the time required for this rill to form was defined as "Trill". Secondary, or tributary, rills formed and contributed to the main rill, but the water tended to continue to erode the initial rill, rather than form a new main rilL Once rilling had begun on the slope, the water began eroding a channel that receded across the crown from the crest of the slope towards the plywood. The rill on the slope went one direction and the channel on the crown went in the opposite direction, but in line with the rill on the embankment slope.
Recorded times for rill initiation to occur on the levee slope for this setup are shown in Table I. For this setup, the levee slope was 0.389m (15.3") in length. The data for this setup yielded unique results. It was the expectation that a high flow rate would erode the slope more rapidly. This relationship may be true for embankments that are already saturated, however, in this setup the plywood divider prevented the model from becoming saturated until overtopping began. So , as water was flowing on the surface, the model was also in the process of becoming saturated. The water at higher flow rates moved more rapidly over the surface of the soil than through the soil, essentially eroding less erodible unsaturated sand and producing larger values of Trill. Because the water at lower flow rates did not move significantly faster over the soil surface compared to flowing through the soil, the model was able to become saturated and more erodible before the rilling process began and yielded smaller values of Trill. In addition to water eroding a channel on the slope, the rill process also involves the location of a specific area on the slope due to geometric or compositional variations, unique to each embankment. To ensure the experiments accurately simulated rill erosion processes, a full-levee model was constructed with a plywood core, so that the geometry of the wood would not be the determinant in the location of rill formation. The results of this setup are shown in Figure 3. Water was supplied to one side of the levee and the water level was allowed to slowly rise until overtopping began. The full-levee model also offered the benefit of more realistic saturation conditions, as the water could flow through the sand to the floodplain side of the levee. Once the water had crossed the levee crown, rill initiation began, at a location influenced by the levee itself, and a channel was eroded on the slope, as well as across the levee crown. Additionally, the use of a full -levee model allowed for a more complete rilling process, as the channel that receded across the levee crown had the opportunity to reach the water-side slope, thereby breaching the levee.
Recorded times for rill initiation and times for the water to cross the 0.203m (8") levee crown, as well as initial (Wi) and final (wr) moisture contents of the levee soil are provided in Table 2. For this setup, each levee slope was 0.259m (10.2") in length.

Full-Levee Setup -Nevada Sand -Kaolin Clay Mix
The full-levee testing of the sand-clay mix followed the same procedure as the full -levee testing with the laboratory sand. However, the preparation of the model varied slightly. The sand-clay levee model did not incorporate a low-permeability core of any kind during the testing. The Nevada sand and kaolin clay also were carefully measured in predetermined proportions based on the sand versus clay content (100-0, 90-10, 85-15) for the particular experiment. Specific volumes of water were measured in order to be mixed thoroughly in the mixer at the desired initial moisture content of 7.5%.
The observed macroscopic erosion processes in the sand-clay mixes were very similar to the erosion process·es observed in the laboratory sand experiments. The water began to saturate the soil while the water level rose on the waterside of the levee, then progressed over the crown of the levee and eventually formed a rill on the lands ide slope and began to recede across the crown, as in the sand experiments. The microscopic erosion processes showed larger clumps of soil (approximately 1.59mm (1 /16") to 3.18mm (I / S ") in diameter) being removed from the levee crown and a type of undercutting taking place, resulting in a faster breach time. The breach time was defined as the time from the beginning of the rill initiation to the time when a channel had completely receded length of the levee's crown. Figure 4 shows the time required for the initial rill ("T_rill") to form on the landside slope of the levee and for full breaching ("T_breach") to occur as a function of kaolin clay content.  As depicted in Figure 4, the faster breach times occur with increasing clay content of the levee and increased size of sediment being eroded from the levee's crown. Conceptually, the data shown in Figure 4 seem counterintuitive as increasingly cohesive levees should require more time for erosion to occur, due to the lower permeability of the material. Physically, the data shown are reasonable if the sediment clumping is considered. Larger and therefore heavier particles were removed from the levee when the rilling process commenced. The flowing water could not carry the mass and deposited the clump of soil farther up on the levee slope, so the clumps of soil accumulated at the leading end of the rill thereby creating a large amount of soil to be eroded down the slope and resulting in longer initial rill formation times. As the channel receded across the levee crown, large clumps of soil were eroded and carried a short distance while suspended in the water before reaching the bottom of the channel when rolling could occur. Because the direction of the erosion occurring on the levee ' s crown is opposite the direction of flow, deposition of sediment did not impede the erosion process, resulting in more rapid breach times as clay content, and clump size, increased Also, because the water travels at a greater velocity on the surface than through the clayey levee, more surface erosion is observed in a shorter period of time. Further testing is required to determine if this trend will continue as the levee becomes dominated by clayey soils.

CONCLUSIONS AND FUTURE CONSIDERATIONS
The model levees eroded more rapidly when fully or nearly fully saturated. A low-permeability core in the center of a levee prevented failure from seepage, and extended the time required for a full breach of the model to occur. A critical clay content in the levees composed of sand-clay mixtures existed at approximately 15 -20% kaolin clay content. A requirement for a core at the center of the levee could be imposed for soils of this composition.
The physical modeling capabilities allow for layered models, such as the inclusion of soil cores, and complex geometries with different crown widths and slope inclinations. Using a geotechnical centrifuge, erosion tests will be performed that will allow simulation and understanding of structures that will be subjected to stresses and forces encountered in earth embankments in the field. Measurement of flow velocity and hydraulic shear stress will be incorporated in future testing. Change detection software will be utilized to gather and process data for multiple layers of soil, allowing for simulation of more complex soil models in the software.