Simulations of air and water flow in a model dike during overflow experiments

Flow in flood dikes, earth dams, and embankments occurs in variably saturated conditions, with pores of the earth material filled partly with water and partly with air. In routine engineering analysis, the influence of pore air is neglected and the air pressure is assumed equal to atmospheric. In some circumstances, for example, during overtopping of the dike by water, the effect of pore air on water flow and stability of the structure can be important. These features cannot be captured with the commonly used Richards equation. In this paper, we analyze earlier experiments on the overtopping of a model dike made of fine sand. During the experiments, a significant amount of air was trapped near the outer slope of the dike, which later escaped through a fracture formed in wet sand. The observations were compared with numerical simulations using the Richards equation and the two-phase immiscible flow model. The deformation and damage of the dike were not modelled, but the initial evolution of the entrapped air pressure (before damage occurred) was in a good agreement with two-phase flow simulations.


Introduction
Water flow in the shallow, partially saturated soils is routinely described by the Richards equation. This equation is based on the assumption that the air (gas phase) in soil pores is connected with atmosphere and can freely escape from the pore space. Therefore, the pore air does not influence the water flow. These assumptions are often used in seepage and stability analysis of earth dams, flood dikes, natural and engineered slopes, and similar earth structures. However, it is well recognized that in many cases, the pore air has a nonnegligible effect on water flow and mechanical behavior of porous material (e.g., [13,14,20,34,35,54,55,62]). During infiltration on large areas, where there is no immediate escape route for the pore air, the advance of the wetting front is slowed down with respect to the solution predicted by the Richards equation [10,16]. Similar effects were observed during infiltration in closed columns [45][46][47]. Moreover, it has been noted that compression of air in pores and the corresponding increase of air pressure may lead to rapid burst movement of air toward the soil surface. This phenomenon is sometimes accompanied by significant damage to the soil structure, with possible dangerous consequence to the soil stability [3,23,61].
In real-life conditions, water infiltration in a hydrophilic porous medium (such as the majority of soils or rocks) nearly always leads to air trapping, i.e., creation of isolated pockets of air, surrounded from all sides with water. Residual air bubbles are often observed at the pore scale, because water preferentially flows as films along the pore walls and portions of air occupying central parts of larger pores can be snapped off [27,28]. Porous media with occluded gas bubbles in liquid-filled pores are often termed quasi-saturated [12]. Gas bubbles in pores can also occur due to cavitation (e.g., [15]). In some situations, air (or other gases) form larger bubbles, with dimensions exceeding several times the average size of soil grains and regular pores-these so-called gassy soils were studied, e.g., by [56,57] and [31]. Trapping of air (or non-wetting fluid in general) is also observed at larger (macroscopic) scale, where the porous medium can be considered a continuum. This can be related to the presence of material heterogeneities, such as layers or coarse inclusions or lenses (having low air-entry pressure) surrounded with finer material (having higher air-entry pressure) [22,32,39,48,52]. In this case, during capillary dominated flow the fine material outside inclusions becomes saturated with water at an elevated value of the capillary pressure, for which the water saturation in coarse material is still low. The air remaining in inclusions cannot escape unless its pressure exceeds the entry pressure of the background material. As a result, local unsaturated zones are formed, with continuous air phase in the pores of coarse inclusion [11,36,40,42].
Macroscopic air trapping can be observed in homogeneous media as well. Such a situation occurs when water infiltrates the domain from all sides, leaving no continuous path for air to escape. This occurs, for example, in flood dikes or earth dams during extreme flood events, when water flows over the top of the structure [43,44]. In order to investigate this phenomenon, a series of experiments were carried out in the Institute of Hydro-Engineering of the Polish Academy of Sciences in Gdańsk. They included infiltration in closed columns [58][59][60] and infiltration and overflowing on a model sand slope [1,5,24], and on a model sand dike [4,23,24]. Overflowing experiments on various model structures in 2D flumes showed formation of air pockets and subsequent damage of the structure, which was attributed to the increase of air pressure.
In this paper, we describe numerical modeling of one of the experiments performed on a model dike. We focus on the phenomenon of air trapping and the associated increase of the air pressure in pores, which ultimately led to a damage to the soil structure. In particular, we attempt to verify if the pressure can be represented with a standard two-phase (air and water) flow model, where parameters are predicted from knowledge about the material only. For the sake of comparison, we also perform simulations with the Richards equation. While the Richards equation is inherently unable to describe changes of the air pressure, it is widely used in studies related to the seepage in dikes (e.g., [19,29,49] and is the basis of many software packages available to the scientific and engineering community, e.g., VS2DTI [18,38], Plaxis [6], or HYDRUS [33]. In the following sections, we describe the experiment (Section 2) and the numerical model (Section 3). Next, we present the method used to obtain material parameters, i.e., the water retention and relative permeability functions of the soil materials (Section 4). Simulation results are shown in Section 5. A sensitivity analysis was carried out in order to assess the influence of different hydraulic parameters. Finally, concluding remarks are given in Section 6.

Experiment
The experiment was carried out in the laboratory of the Institute of Hydro-Engineering of the Polish Academy of Sciences in Gdańsk. The main part of the experimental setup was a Plexiglas flume having dimensions 200 × 100× 4.5 cm (Fig. 1) which allowed to create twodimensional fluid flow and two-dimensional plane strain conditions for soil deformation.
Two types of fine sand were used in the physical model, FSa1 for the bottom layer and FSa2 for the main structure. Both types of sands were well sorted, with steep granulometric curves (Fig. 2). The structure was created by  the sandy rain method, i.e., gently pouring the sand without compaction, using a specially designed device. Consequently, the obtained angle of the slope was very close to the friction angle of the soil, which was equal to 34 • for FSa2 and 37 • FSa1. The porosities were respectively 0.55 and 0.52 [4,25].
As the first step of the experiment, the bottom layer of sand was saturated and water table established at the base of the dike (Fig. 1). After that, an overtopping simulation was performed. The water table at the left-hand side of the dike was raised to the level of the crest within a period of 480 s, allowing water to flow over the top of the dike. In order to stabilize the dike, which was undergoing significant displacements, the water table was raised also at the righthand slope until its level reached half of the height of the dam (this process took 60 s).
As a consequence of overflowing, a large unsaturated area surrounded from all sides with saturated or quasisaturated material occurred in the inner part of the dike, near the right-hand slope. The air pressure in the entrapped air zone was measured by a pressure sensor. Its position was chosen based on preliminary overflowing tests (Fig. 1). The measured air pressure values (relative to the atmospheric pressure) are shown in Fig. 3, with pictures showing the saturation distribution in dike at the different stages of the experiment.

Fig. 3
Pressure change in unsaturated area: 1, before overtopping; 2, just after overtopping; 3, state corresponding to the max. value of pressure in the area; 4 and 5, decreasing of the air pressure caused by the increasing volume of the macropore [23][24][25] A slight pressure buildup in the air phase was observed since the beginning of the infiltration. In the eighth minute of the test, water flowed over the crest of the dike and air pressure raised up as the connection between closed unsaturated pore space and atmosphere was lost because of a thin layer of flowing water on the slope and a layer of saturated sand at the surface of the slope. Due to the imposed boundary constraints (Fig. 1), the left-hand slope and lower part of the right-hand slope submerged with water), the overall pattern of the pore air flow was toward the upper-right part of the dike. Further infiltration of water caused reducing the volume of the unsaturated area and noticeable increase in pressure (Fig. 3, between stages 1 and 2). A void space (macropore) appeared at the boundary between the unsaturated and (quasi-)saturated zones in the form of a cavity surrounded by porous material and occupying the whole thickness of the model. Its volume increased as the air entrapment zone was reduced. Then, there was a drop in the pressure due to air escaping through the channels (cracks) appearing in the thin saturated layer in the upper part of the right-hand slope. Subsequently, due to ongoing infiltration and deformation of the slope surface, porous channels were closed, and as a result, the air pressure rose again. The maximum was reached in the eleventh minute when a narrow crack appeared again, linking the unsaturated area to the outer slope. The volume of this crack grew with a simultaneous drop of the air pressure (Fig. 3, between stages 4 and 5).
The cracks appeared when the air pressure measured by the sensor is in the range of 400 to 650 Pa. This corresponds to the weight of a layer of 2-3 cm of water-saturated sand (assuming density of about 2 g/cm 3 ). While this is a very crude approximation, it shows that the increased air pressure might be sufficient to reduce the effective stress in the soil skeleton near the surface of the left-hand slope and cause the damage to the soil structure.

Numerical model
Numerical simulations of the experiments were carried out using the two-phase immiscible flow model and the Richards equation. Since our models do not take into account the deformation and damage of the soil skeleton, we focus on the part of experiments before the appearance of cracks. We assumed that fluids are immiscible, soil particles are incompressible, the process is isothermal, the soil skeleton is rigid, and has isotropic permeability. Under these assumptions, the governing equation can be written in the following form (e.g., [17]: where α means air (a) or water (w) phase, respectively, ρ α is the density (kg/m 3 ), S α is the saturation (S a + S w = 1) (-), n is the porosity (-), t is the time (s), κ is the intrinsic permeability (m 2 ), μ α is the dynamic viscosity (Pa s), k rα is the relative permeability (-), p α is the pressure in fluid (Pa), g is the gravity (m/s 2 ), and z is the depth (m). In the twophase flow model (2PH), we have two separate equations, where in the first one index α is replaced by "w" and in the second one by "a". The flow model proposed by Richards (RE) is obtained by retaining only one equation, where the considered phase is water. Both fluids are assumed to be compressible, with water compressibility coefficient of 4.5×10 −10 Pa −1 . The relationship between density and pressure of air is taken into account by the ideal gas law. The density of pore air varies according to the changes in the content of water vapor, but these effects are neglected in the present work. The pressure difference between the non-wetting fluid (air) and the wetting fluid (water) is called capillary pressure or suction. The value of the suction pressure is related to the water saturation. This relationship is called soil water characteristic curve (SWRC) or soil water retention curve (SWRC) and is required as a model input together with functions describing the dependence of relative permeabilities on saturation. Here, we use the well-known formulations of van Genuchten-Mualem [30,50]: and Brooks-Corey-Burdine [7,8]: where S ar is the residual saturation of fluid (-), S ew is the effective water saturation (-), p g is the scaling pressure (inverse of the parameter α of the original van Genuchten model) (Pa), n g (-) and m g (-) are the exponents in the retention function, p e is the air entry pressure (Pa), λ is the shape parameter (-). In contrast with the van Genuchten model, the Brooks-Corey capillary curve includes a distinct value of the air entry pressure. The governing equations for unsaturated and two-phase flow were solved using an in-house code implemented in Fortran. The algorithm is based on the control-volume finite-element approach, which combines features of the finite volume method (FVM) and the finite element method (FEM)-for details, see, e.g., [17] and [37]. In order to avoid the differences between solutions, all numerical simulations described in this paper were carried out using the same approach to discrete 2PH and RE mathematical models. The algorithm is based on vertex-centered finite volume discretization in space (Fig. 4). In this case, we used quadrilateral elements as the primary grid, on which the dual grid of finite volumes was superimposed. The gradients of water and air pressures through each segment of the control volume boundary are approximated using shape functions for the corresponding finite element and the relative permeabilities were calculated as upstream values. All simulations were performed on the same numerical grid consisting of 779 nodes. As shown in Fig. 4, the grid included the upper structure of the dike and a major part of the base (we assumed that flow in the outermost regions of the base can be neglected).
For time discretization, we used the fully implicit firstorder finite difference method. In each time step, the discretization results in a system of nonlinear algebraic equations, which must be solved with respect to the primary unknowns. In the case of two-phase flow, we choose the water pressure and water saturation as the primary variables, in accordance with the most common pressure-saturation formulation. For RE, we have only one unknown per node, which is the water pressure. The nonlinear algebraic system arising at each time step are solved using Newton-Raphson iterative scheme with line search. The time step varied in the range from 6.2×10 −4 to 3.147 s. It was adjusted automatically based on the performance of the iterative solver.
The boundary conditions varied in time according to Fig. 1. At each stage of the experiment, we imposed hydrostatic value of the water pressure p w and S w = 1 on the submerged parts of the slope, while on the dry parts we assumed no water flow and atmospheric air pressure. The base and vertical sides of the flume were considered impermeable for both water and air. Initially, the area above the water table is in unsaturated state. Water pressure in the domain is assumed to vary hydrostatically from 5.39 kPa on the bottom to − 4.14 kPa on the top. Overtopping is represented by a time-dependent pressure condition. For example, for the top boundary in the beginning of the simulation, we assume that the water pressure is negative (− 5395.5 Pa which is equal to the product of water density, gravity, and dike high) but after overtopping this value changes to 0.

Estimation of material parameters
Since no measurements of the SWRC or permeability for the sand materials were available, we estimated their hydraulic parameters from the granulometric curves. Such an approach is well established in the literature, although it cannot be considered very accurate. A number of formulas, known as pedotransfer functions, were proposed to compute hydraulic function parameters from basic soil properties, such as granulometric curve, bulk density or porosity. Here, we used the modified Kovacs (MK) method, developed by [2].
The MK model assumes that water is held by capillary forces, responsible for capillary saturation (S wc ) and adhesive forces, causing saturation by adhesion (S wa ). The proposed equation is expressed as: where denotes the Mackauley brackets x = 0.5(x + |x|). S wc and S wa depend on ψ (suction expressed in centimeters of H 2 O): S wa = a c C ψ (h co ) 2 / 3 e 1 / 3 (ψ) 1 / 6 (10) where h c (cm H 2 O) is the equivalent capillary rise, described below, m and a c are dimensionless parameters (for granular and soils m = 0.091 and a c = 0.013), ψ 0 is the maximum suction in oven-dry conditions, equal to 10 7 cm H 2 O and ψ r is related to h c :  The equivalent capillary rise h c is defined as: where σ w = 0.073 N/m is the surface tension of water, β w = 0 is the contact angle between water and solid phase, α = 10 is a shape factor, e is the void ratio, and D H is an equivalent particle diameter, which can be approximated using the following function [2]: where d 10 is diameter corresponding to 10% passing on the cumulative grain size distribution curve and C U is the coefficient of uniformity (C U = d 60 /d 10 ).
The MK retention functions of both sands, obtained in tabular form, were converted to VGM and BC functions using the RETC computer program [51]. The residual saturations were assumed S rw = 0 for water and S ra = 0.02 for air. The fitted parameters p g and n g are given in Table 1, and the curves are shown in Fig. 5. The permeability coefficient can also be estimated based on the granulometric curve. There is a large number of available formulas, developed to estimate permeability with respect to water in saturated conditions (k ws = κρ w g/μ w ), from which the intrinsic permeability κ can be also calculated. For a review of such methods, see, e.g., [9,41,53]. Here, we used four well-established formulas. Kozeny-Carman and Zunker equations were taken from [41]. They include dependence of k ws on porosity. They also require a representative grain diameter, which is calculated as an average value based on the whole granulometric curve. The Hazen and USBR equations are respectively k ws = 1200 · d 2 10 and k ws = 311 d 2.3 20 . where d 10 and d 20 are diameters corresponding to 10 and 20% content of grains in Fig. 2.
The calculated values are listed in Table 2.
The differences in computed permeability values are significant, but such discrepancies between empirical formulas occur commonly (e.g., [41]). For the core of the dike (FSa2), the largest value was obtained with the Hazen formula and the smallest values with the USBR formula. After 540 s, significant differences are observed between the experimental results and computations. This is due to the escape of air through macropores, which leads to oscillations and reduction of the measured pressures. This effect cannot be represented in our numerical model. Note that oscillatory behavior of the pore air pressure and related air flux between subsurface and the atmosphere was modelled by [10] using the leakage concept, where an additional set of parameters was assigned to the ground-atmosphere interface. Such an attempt was not undertaken here, and it may be a topic of further study.
Numerical solutions start to diverge significantly approximately at the same time as the air breakthrough occurs in the experiment. This time corresponds to the full development   of an entrapped air region inside the dike, separated from the right-hand slope by a layer of quasi-saturated sand. In the numerical model, further infiltration causes air compression, proportional to the permeability of the medium. Note however, that even though the numerical simulations do not reproduce the measured pressure drop, there is a visible change of shape in the pressure evolution curve at the breakthrough time. This corresponds to the start of air outflow through the quasi water-saturated layer, which developed along the surface of the slope during overflow. For the largest permeability estimate (the Hazen formula), a good agreement with measurements is obtained again for the second peak value of the air pressure. However, the numerical simulation predicts further buildup of the air pressure up to the final value of 1605 Pa (results not shown here). However, these high values of the air pressure correspond to the state of residual air saturation and are determined by the elevated values of the water pressure, imposed by the external water table. Other estimates of permeability lead to significantly lower rate of the pressure increase. Figure 8 shows the distribution of water saturation, pore air pressure, and normalized pore air pressure at time 580 s, shortly after the appearance of cracks. The air pressure is calculated from the primary variables of the 2PH model, p a = p w + p c (S w ). The largest air pressure values occur in the lower-left part of the dike. In this region, air is at residual saturation, with capillary pressure close to 0, but the water pressure is high, because of the position of the water table, which means that the calculated air pressure is also high. We also show the "normalized" air pressure, i.e., the air pressure multiplied by the effective air saturation S ea = (S a − S ra ) /( 1 − S ra ). According to many authors [21,26], this corresponds to the part of the effective stress in soil attributed to the air phase. The highest values of the normalized air pressure occur in the region of entrapped air in the upper part of the right-hand slope.
Sensitivity studies were also performed with respect to the parameters of the water retention curve, i.e., van Genuchten parameters n g (Fig. 9) and p g (Fig. 10) and Brooks-Corey parameters n b (Fig. 11) and p e (Fig. 12). It can be seen that the pressure scaling parameters have much larger influence on the evolution of the air pressure than the shape-related parameters. Increase of the p g or the p e allows for development of larger air pressures in the entrapped region. In such a case, larger difference between air and water pressure is necessary, in order that the air could flow through the quasi-saturated layer along the slope. While the discrete air entry pressure is included only in the Brooks-Corey model, a similar effect is observed for the p g of the van Genuchten model. In both cases, changing the initial parameter value by about 12 to 15% caused very clear changes in the solution. It can be seen also in Fig. 13 where the saturation maps for different values of scaling pressure and air entry pressure has been compared. In contrast, the parameters n b and n g , related to the shape of the water retention curves, did not have significant effect on the results, even though their relative changes were larger (15 to 20%).
For the purpose of comparison, we also performed a numerical simulation using the Richards equation and the Brooks-Corey-Burdine parameters reported in Table 1. Figure 14 shows the water saturation distributions in the dike structure at time t = 660 s, corresponding to the maximum air pressure measured during the experiment. The presented saturation fields were calculated with the Richards equation and the two-phase model (using the same parameters). Note that due to deformation during overtopping the model dike slightly changed its dimensions (it became a little lower and wider). There is a significant difference in the extent of the unsaturated area which is smaller in the RE solution than in the 2PH. However, as discussed above, the saturation distribution is affected by the choice of the pressure scaling parameters. It is possible to obtain a good match with the experimental results using RE with air entry pressure reduced by approximately 15% or hydraulic conductivity reduced by approximately 25%. Such variations of parameter values are within the expected range of uncertainty. Thus, the evolution of saturation could be represented either by the RE or 2PH models, but obviously only 2PH model is capable of reproducing the increase of air pressure in pores.
The results from 2PH model are in good visual agreement with the experiment, although it can be noted that the sharp interface between the saturated and unsaturated area, which formed parallel to the right-hand slope of the dike after the start of overflowing, is not correctly represented. It seems that in the model, the pore air can escape more easily from the dike than in the experiment. One possible reason might be related to neglecting hysteresis in the water retention function. The upper layer of the slope, which is saturated by water after overflowing becomes later invaded by the air escaping from the inner part of the dike. So, hysteresis might play some role, especially in the surface layer of the slope, leading to an increased value of the air entry pressure. On the other hand, the position of the zone of air outflow on the dike slope is correctly reproduced by the model. In the numerical simulations, the air escape is modelled as a continuous process, while in the experiments, discrete outbursts were observed, connected with damaging soil structure and creation of macropores.

Conclusions
The research described in this paper provided insights on the possible role of the air phase during overflowing of dikes and similar structures. During the overflow experiment on a model dike, a large amount of air was trapped which finally led to the formation of macropores and cracks. The general pattern of water infiltration observed in the experiment was in qualitative agreement with earlier numerical simulations for real-sized dikes [43,44]. In this work, an attempt was undertaken to reproduce the experimental results with numerical simulations based on a two-phase flow model with parameters estimated from the granulometric curve using the modified Kovacs method by [2]. The agreement with the experiment was good for the initial stage of infiltration, especially when the Brooks-Corey-Burdine model was used together with the intrinsic permeability calculated according to the Hazen formula. Due to the limitations of the applied model, it was not possible to reproduce the later oscillations and dissipation of the pore air pressure, related to the formation of macropores. In this regard, further studies are necessary. Future modeling efforts could include coupling between flow and deformation model, introduction of the air escape mechanism similar to the one described by [10] and/or hysteresis of the water retention curve.
Sensitivity analysis showed that in the earlier stage of infiltration, before the formation of macropores the evolution of the air pressure is not sensitive to the changes of the permeability and the retention function parameters. In the later stage, the air pressure values are sensitive to the pressure scaling parameters of the Brooks-Corey and van Genuchten retention curves and to the intrinsic permeability value. The shape parameters n g and n b influence the results to a much lower degree The changes in pore air pressure seem to have an important influence on the mechanical behavior of the structure. Their evolution cannot be reproduced with the Richards equation. Also, in terms of the distribution of the water saturation, the two-phase model was closer to the experimental results than the Richards model. However, it was possible to reproduce the saturation pattern by changing the hydraulic parameters in the Richards equation. These results highlight the importance of using multiphase flow models for mechanical analysis of structures made of earth materials, confirming observations from other authors (e.g., [20,34,35,62]).