Numerical simulations of overflow experiments on a model dike

. Dike failure due to overtopping is one of the important factors, which should be considered in the dike designing process. Although the overflow is characterized by the relatively low risk of occurrence, in many cases dikes are totally destroyed or seriously damaged. An interesting phenomenon occurring during overflow is the trapping of air in pores of the unsaturated soil material. As the infiltration progresses from all sides, the air pressure in the unsaturated region increases, which may ultimately lead to damage of the dike structure. It happens when the air is expulsed in form of bursts and forms large macropores. Such a behaviour evidenced in laboratory experiments. In this study we attempt to simulate the evolution of stress field in the model dike subjected to overtopping. The results are in qualitative agreement with observations, showing that formation of the first macropores occurs in the direction perpendicular to the minor principal stress in the soil mass along the dike slope edge.


Introduction
Dikes are usually build along the rivers to protect inhabited areas from floods. The existence of the dike retaining water is always connected with the risk of their failure which can be disastrous and cause huge economic losses. So in the design process it is necessary to take into account various types of failure mechanisms. Unfortunately assumption of the full protection for this type of structure is almost impossible because there is always a risk, even if a very small one. It is associated with the occurrence of phenomena such as overflowing, internal erosion or sabotage [1].
In general, the causes of dike failures can be divided into three groups: the loss of stability (shearing failure along the sliding surface or cracking due to settlement), erosion, damage caused by water and overflowing. The statistics published in [1] showed that occurrence of the overflow caused approximately 40% of all dike failures. A number of authors have studied the mechanism of dike failures due to overflowing [2,3,4,5,6,7,8]. They are in agreement that fast raise of the water table level cause air trapping on pore and structural level. In the first case small bubbles of air are more or less uniformly dispersed in water-filled pores of the material. In the second case a large unsaturated (i.e. mostly air-filled) area surrounded by the water-saturated material is created. Both phenomena are connected with the loss of macroscopic homogeneity which is particularly undesirable in dike structures [8,9,10,11]. This problem is difficult to investigate in-situ due to the lack of appropriate monitoring of the water saturation. Therefore an important direction of research are experiments and numerical methods which allow to simulate dike behaviour in controlled conditions. In this work we analyze the results of laboratory experiment performed at the Institute of Hydro-Engineering of the Polish Academy of Sciences in Gdańsk. The experiment was carried out on a two-dimensional model dike in small scale, in order to observe the air trapping during the overtopping process. The details of the experiments are provided below. One of the interesting findings was the appearance of large voids (macropores) during air entrapment ( Fig. 1 and Fig. 2), resulting from an increase of the pore air pressure in the unsaturated zone (up to about 0.65 kPa over the atmospheric pressure). This type of void was also observed in earlier experiments, carried out on a slightly different dike configuration [6,12]. Numerical simulations [12] showed that the evolution of air pressure in the model dike can be successfully represented with a two-phase flow model (accounting for both air and water phases). However, the stress state of soil was not considered in the simulations, while it clearly had an important influence on the formation of cracks. It is generally assumed that cracks induced by increasing pore pressure occur in the direction normal to the minor principal stress [13]. The in-house code for two-phase flow used in [13] does not account for stress and deformation analysis. In order to gain at least qualitative insights into the stress field evolution, in this paper we perform numerical simulations of the overflowing process with ZSoil code [14]. ZSoil is capable of coupled analysis of soil deformation and unsaturated water flow, although the air flow is not modelled. Despite this limitation, the obtained information on the total and effective stress field can be useful for the analysis of the macropore formation.

Experiment
The experiment was carried out at the Institute of Hydro-Engineering of the Polish Academy of Sciences. The model dike was set up in a Plexiglas flume having dimensions 200x100x4.5 cm which allowed to create two dimensional flow and deformation conditions (plane strain). During experiments no signs of preferential flow of flume walls were observed. Two types of sands FSa1 (bottom) and FSa2 (main structure) were used to form the model dike. The sand was poured down from a special device, without further compaction, so the slope of the dike was close to the natural angle of repose. As the first step of the experiment the bottom layer of sand was saturated. After that overtopping was introduced. Water table at the left side of the dike was raised to the crest in 480s which caused the water overflow on the right slope. This stage has been hold by 60s with simultaneous rise of the water table on the right side of the dike. The mentioned procedure was necessary to prevent collapse of the structure. The last step was lowering of the water table on the left side to the end of the simulation which occurred at 2000s. The geometry and boundary conditions of the experiment are presented in Fig. 4.
During overflowing an unsaturated zone of entrapped air appeared in the dike structure (Fig. 1). The air pressure in the entrapped zone was measured by a pressure sensor, shown in Fig. 4. Increase of the pore air pressure led to the damage of soil structure and creation of macropores. The first type of void appeared between the saturated and unsaturated zones in the form of a cavity. Its volume increased as the unsaturated soil region was reduced and simultaneously the air pressure in unsaturated sand increased. After reaching the pressure value of about 400 Pa (over the atmospheric pressure), the gas phase escaped, creating another macro-pore, which linked the existing void with the slope of the dike. This was accompanied by a drop in the air pressure. Later the air pressure rose again to the value of about 0.65 kPa and finally decreased to the atmospheric value as the whole soil in dike structure became water-saturated [12].
The basic soil parameters such as friction angle, density and porosity where taken from [6], hydraulic parameters describing water retention function were determined using RETC program [18] (tab. 1) while filtration coefficient was calculated based on analytical formula proposed by Hazen: 2 10 1200 Elastic Young's-module was taken as 60 MPa for both types of sands. This value was considered representative for the medium dense sand in the model dike structure (in contrast to the real dikes where the soil is heavily compacted). The Poisson's ratio was calculated using the following equation based on effective friction angle: where w is phase occupying the pore space (here water), w -density, Swsaturation, nporosity, t -time, intrinsic permeability, μw -dynamic viscosity, krwrelative permeability, pwpressure in fluid, g -gravity, z -depth and Cwpstorage coefficient: where w is volumetric content of water, w is water compressibility coefficient. Equation (3) contains several dependent variables, so it has to be completed by additional constitutive relationships. Capillary pressure is equal to the negative water pressure. The value of the capillary pressure is related to the water saturation. This relationship is called SWRC (soil water characteristic curve) or SWRC (soil water retention curve) and can be expressed by various analytical formulae which examples of which can be found in [15,16,17]. In ZSoil the SWRC is implemented as a simplification of the van Genuchten formula with fixed values of n and m exponents, equal to 2 and 0.5 respectively: where Swr is the residual water saturation, Sw -water saturation, pg -scaling pressure in the retention function and pc -capillary pressure. Permeability of the medium is obtained by scaling the tensor for fully saturated medium by a scalar value function dependent on the saturation ratio: Due to the coupling of the water flow with deformation it is necessary to extend the effective stress concept to account the partially saturated medium. In ZSoil code a following formula has been adopted: where ' -effective stress,  -total stress. Water saturation equal to 1 denotes fully saturated state, between 0 and 1 is partially saturated state and 0 corresponds to dry conditions.

Discretization and boundary conditions
ZSoil is written in the MS-Windows environment and it is based on the finite element method with linear, quadratic or cubic shape functions for spatial discretization and fully implicit backward finite difference formula for the integration in time. Evaluating of all integrals resulting from the problem formulation it is necessary to use numerical integration. In ZSoil a Gaussian quadrature is used. The nonlinear algebraic systems arising at each time step can be solved using Newton-Raphson iterative scheme with backtracking or BFGS algorithm. It is possible to use steady state or transient analysis of the water flow for partially or fully saturated medium with time dependent boundary conditions. In the case of single phase water flow, the fluid pressure is a basic unknown. When considering the two-phase model additional unknown is nodal displacement. An important advantage of the software is that user can change the parameters of the iterative solver, for example maximum number of iterations per step. To describe soil deformation an elastic model has been used in all calculations presented in the article.
All simulations were performed on the same numerical grid consisting 951 nodes and 870 quadrilateral elements. As shown in fig. 5, the shape of the dike base was assumed different than in experiment because it speeded up the calculations while not affecting the results. The numerical model was implemented in four stages. The first three represented respectively: placement of the base, placement of the structure of the dike and saturation of the base. The last stage simulated the overflowing with boundary conditions varying according to fig. 4. At each time step we imposed hydrostatic value of the water pressure pw on the submerged part of the slopes, why on the dry parts we assumed no-flow conditions. The bottom edge of the dike was considered impermeable for water and fixed with respect to displacements in both directions. Before rising the water table (end of the third step) the whole area above the water table was in unsaturated state. Water pressure in the domain was varying hydrostatically from 1.10 kPa at the bottom to -5.50 kPa at the top.  Fig. 6 shows that the regions of the minimum values of the effective stress (in horizontal direction) are located inside the dike near the upper part of the right-hand slope. This is the region where formation of cracks occurs (compare fig. 1). We expect the actual values of the effective stress to be different from the computed ones, since the numerical model did not take into account the influence of the pore air pressure. However, one should be aware that the definition of effective stress in a porous medium saturated by multiple fluids is not straightforward and is a subject of ongoing research [18]. Fig. 7 presents the evolution of principal stress components in the considered region of the dike. It can be seen that the ratio of the minimum to maximum effective stress diminishes as the infiltration proceeds, while the direction of the minor principal stress is approximately normal to the right-hand slope of the dike. This is in agreement with the observed macropore formation, since the void space has its larger dimension oriented approximately parallel to the slope, i.e. parallel to the major principal stress and perpendicular to the minor stress. This is consistent with the theoretical analysis of fracturing due to pore pressure increase [13].  The calculated values of displacements and water pressure at a node located in the region of crack formation are shown in figs. 8 and 9. Since the numerical model was implemented in four stages, the plots present only results from the fourth stage (overflowing). Fig. 8 shows the evolution of the water pressure for the simulations with and without soil deformation. Until 550s (end of rising water table on the right side and beginning of decreasing the water table on the left side) both curves are in agreement. When the process of rapid infiltration from both sides ends, the differences appear. They can be explained by the fact that including soil deformation of the model adds additional retention capacity to the soil medium, which slows down the propagation of the pressure changes. The displacements in both directions ( fig. 9) are relatively small which is inconsistent with experiment where significant deformation and damage of soil structure has been observed. This is because we assumed a simple elastic model for our calculation, which does not account for plasticity or brittle failure.

Conclusion
Despite the limitations of the numerical model used in this study, the calculated evolution of the effective stress field is in the qualitative agreement with the soil behaviour observed in the experiments. In the region where soil cracking occurred primarily the values of the effective stresses are the lowest. Additionally the minor principal stress is much smaller than the major stress. The minor stress is oriented perpendicular to the slope and approximately perpendicular to the longer dimension of the macropore observed in the experiment. This is in agreement with the theory of fracturing as discussed by [13]. Improved analysis of this phenomenon, based on the more sophisticated flow and deformation models, is warranted and will be the subject of further work of the authors.