Influence of the air phase on water flow in dikes

Numerical models are often used to describe flow and deformation processes occurring in dikes during flood events. Modeling of such phenomena is a challenging task, due to the complexity of the system, consisting of three material phases: soil skeleton, pore water and pore air. Additional difficulties are transient loading caused by variable in time water levels, heterogeneity of the soil or air trapping. This paper presents a brief review of the influence of the air phase in soil on water flow and pore pressure generation, with focus on applications related to stability of dikes, earth dams and similar structures. Numerical simulations are carried out to investigate the differences between the Richards equation and the two-phase flow model, using an in-house code based on the finite volume method. A variety of boundary problems are considered, including seepage through flood dikes, dike overtopping and water level fluctuations. Special attention is paid to the problem of air trapping, which occurs when water flows over the top of a dike. Such a phenomenon occurred during experiments on model dikes reported in the literature, ultimately leading to development of cracks and damages in dike structure.


Introduction
When considering water movement in soil, it is often assumed that the air present in the soil is connected with the atmospheric air, which has approximately constant pressure.It means that the air in pores can escape freely before advancing water and thus has no influence on the infiltration rate.In most flow scenarios considered in civil and environmental engineering the flow of the air phase is neglected.Consequently, most of the available software for modeling unsaturated zone flow is based on the Richards equation, which does not account for the air flow.However, there is an increasing evidence that in some situations air flow must be taken into account: heterogeneous soils [1][2][3], rapid downward infiltration [4], water table fluctuations [5], seismic events [6] or overtopping of dikes [7][8][9][10].In such cases the two phase flow model should be taken into account.In this type of analysis two equations (air and water) are coupled together.In general, each of the phases consists of multiple chemical components, which can move between phases.Pore air, for instance, is a mixture of gases, including water vapor, while pore water contains many dissolved substances, including gases.The number of phases and components included in the mathematical model depends on the problem under consideration.In many applications focusing on the water flow, a sufficient accuracy can be achieved with a simplified model, where both air and water are considered as immiscible and the deformation of the solid skeleton is treated in a simplified manner [2,11].The main goal of the paper is to present numerical simulations of water and air flow in dikes carried out using an in-house code developed by the authors, based on the control volume finite element approach (CVFE) [2,11].We compare solutions obtained using the full two phase model with the solutions of the Richards equation, which neglects the air flow.

Equations
Mathematical model used in this article is based on the following assumptions: fluids are immiscible, soil particles are incompressible, process is isothermal and all deformations are caused only by changes in the pore pressure and can be described by a linear isotropic elastic law.Under these assumptions, the governing equation can be written in the following form: where D means air or water phase respectively, U D -density, S D -saturation, n -porosity, t -time, N -intrinsic permeability, μ D -dynamic viscosity, k rD -relative permeability, p D -pressure in fluid, g -gravity, z -depth.In the two phase flow model, we have two separate equations, where in the first one index D is replaced by 'w' (water) and in the second one by 'a' (air).The flow model proposed by Richards is obtained by retaining only one equation, where the considered phase is water.Although one of equation is omitted, equation (1) retains its nonlinearity.However it is easier to solve than the two phase flow model.The governing equations contain several dependent variables, so they have to be completed by additional constitutive relationships.The values of the two saturations sum up to one: Pressure difference between the non-wetting fluid (air) and the wetting fluid (water) is called capillary pressure or suction: The density of each fluid is assumed to depend only on the fluid pressure: The value of 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.The most popular models are proposed by Brooks-Corey, van Genuchten or Gardner.In this article we used van Genuchten formula: where S Dr is residual saturation of fluid, S ew -effective water saturation, p g -scaling pressure in the retention function, n g and m g are the exponents in the retention function.
Permeability of the medium describing the Mualem-van Genuchten model: Volumetric deformations of the soil skeleton are accounted for in a simplified way, by assuming linear relationship between the porosity and the average pore pressure: where n 0 is the initial porosity and β s is compressibility coefficient, related to soil Young modulus E and Poisson ratio v : p' is the average pore pressure and can be defined as: Nonlinearity of main equation or equations results from the relationship between saturation and capillary pressure or permeability and type of the soil (clay, sand, silt).Besides definition of the constitutive relationship it is also necessary to select appropriate par of unknowns.We distinguish three most common solutions: pressure formulation, pressure-saturation formulation and saturation formulation.In our previous study on this subject [10] we used a commercial code FlexPDE, where pressure-based formulation was implemented.However, the mixed formulation (pressure-saturation) seems to be better suited for the cases where one of the phases may disappear (such as the air phase during variably saturated flow in soils) [2,11].Therefore the mixed formulation was used in this paper.

Numerical implementation
The governing equations for unsaturated and two-phase flow were solved using inhouse code implemented in Fortran.The algorithm is based on the control-volume finiteelement approach, which combines features of the finite volume method (FVM) and the finite element method (FEM) -for details see e.g.[2].Spatial discretization of the governing equations is performed on a dual grid, where finite volumes (control volumes) are built around nodes of the primary finite element grid.In two dimensions either triangular or quadrilateral elements can be used, while the resulting control volumes are irregular polygons, see Fig. 1.The fluxes of water and air through each segment of the control volume boundary are approximated using shape functions for the corresponding finite element.We use a vertex-centered scheme, where the primary unknowns (water pressure and water saturation) are associated with nodes of the primary finite element grid, which become centers of the finite volume of the dual grid.In such a way the values of unknowns at the outer boundaries of the domain are obtained directly from the solution.For the discretization in time the fully implicit first-order method was used.The nonlinear algebraic system arising at each time step are solved using Newton-Raphson iterative scheme.One of the main advantages of using an in-house code is a greater flexibility with respect to the choice of boundary conditions and material constitutive models, as compared to many commercial software packages.In unsaturated flow simulations it is often necessary to switch between Neumann and Dirichlet boundary condition, for example when modelling the infiltration or evaporation processes at soil surface, or seepage faces on slopes, such as the one presented in the example below.In this article this solution was used in modelling the seepage face of the dam.When water flow out only from some part of the boundary we assume a Dirichlet but for the rest which is unsaturated we assume Neumann.Position of the saturated-unsaturated interface is not known a priori and must be obtained iteratively during simulation.

Geometry, material parameters and boundary conditions
The influence of air phase on the water flow in dikes is investigated for two different patterns of seepage.The first test represents seepage through an embankment caused by an increase of the water level at the outer slope (Fig. 2).We consider homogeneous dike placed on an impermeable base.Soil parameters correspond to loamy sand are listed in Table 1 [13].Initially the dike is in unsaturated state with water pressure equal to 0 (atmospheric pressure) at the base and decreasing hydrostatically to -29,43 kPa at the top of the dike.In the whole domain air pressure was assumed constant and equal to 0 (atmospheric pressure).The time of rising of the water table is assumed to be negligibly small, so the water table level was given at the height of 2,5 m right at the beginning of the calculations.At the submerged part of the slope the boundary condition is given as hydrostatic and rising hydrostatically from 0 at the water table to 24,53 kPa at the bottom of the dam.For the rest of the left slope and the crest is assumed that gas phase can freely escape.The liquid phase can only flow out from the right sight of the dam.This part of the boundary is considered as seepage face.It means that it is impermeable for water as long as the water pressure at the boundary is negative.If the water pressure reaches 0, water starts to flow out from the soil, but no pressure built up is allow due to surface runoff.Boundary condition changes to pressure-type, p w =0.For two phase flow model simulations, a constant value of air pressure p a =0 is assumed along the slopes and top of the dike.
Numerical grid was performed using quadrilateral elements with computational points in the grid nodes (vertex centred method Fig. 1).Here we show results for an 800 elements and 626 nodes because for results for denser grid were qualitatively the same.Time step varied from 10 -12 to 1000 s.

Results
Results of the simulations are shown in fig.3-4.Fig. 3 show the distribution of water saturation in dike at an early stage of the simulation (2·10 5 s).Shape of the saturation front is qualitatively similar but it is definitely longer compared to the 2PH model.At the end of the simulation this difference is a little smaller (steady state has not been reached).Figre 4 shows the total outflow of water from the right slope of the dike (seepage face) as a function of time.Initially value was nearly the same for both models, greater differences can be observed in the middle of the simulated period.It is associated with reaching the seepage face by the water.From this time differences are more visible.This is due to later start of the outflow for two phase model and consequently about 10% smaller discharges at the end.Received results are consistent with expected behaviour of both solutions.For Richards we expected that water should propagate faster than for 2PH.It is caused by the lack of presence the air in the domain, so water can flow out without any difficulties.While in 2PH water needs to displaced air during the whole process therefore infiltration takes a little longer.

Overtopping of a dike
Second simulation shows overtopping of a dike (fig.2).This situation occurs when water table rise over the crest of the dam eg.during heavy rainfall.Geometry, material parameters and initial conditions are the same as in example 1.However, this time we assumed that water rise during 5 hours to the top of the dike and remained there till the end of the simulation.Initially crest and inner slope were treated as impermeable for water.When the overtopping started, the boundary condition was changed to p w =0.As in the first example air pressure is p a =0 for all of the boundaries with the exception of the bottom.Numerical grid and range of time steps were the same as in example 1.

Results
Figures 5-6 show distribution of water saturation shortly before achieving fully saturation according to Richards equation.We can observed that in 2PH model at the bottom right corner is still large unsaturated or partially saturated area while in Richards the dike is almost fully saturated.In the two-phase solution inflow of water is hampered by the necessity to displace the air through the water saturated medium surrounding the unsaturated zone.In Richards equation water infiltrates with the same velocity into an unsaturated area, regardless the possibility of air escaping, thus region of entrapped air is definitely larger in 2PH.Process of saturation during time can be also see on figure 5. We can observe that the fully saturated state occurs about 20 hours later in 2 phase flow model than in Richards.

Conclusions
Numerical simulation presented in the article showed differences between the unsaturated flow model (Richards equation) and the two-phase model.Infiltration process occurring in dike is slower because water needs to displace air which is typically not accounted for in geotechnical and hydrological modelling software [14].In addition when possibility of outflow for air phase is further limited, the impact is even greater.This can happen when we are dealing with overtopping of dike or heavy rainfall during infiltration [15][16][17].The results confirm the findings from our previous study [10] performed using a commercial code FlexPDE based on the standard finite element formulation.The code developed by the authors, described in this paper, offers greater flexibility with respect to the boundary conditions and material constitutive models and therefore is suitable for further extensions.
Problem of air trapping in large soil structures has a significant impact on their stability and surface structure.In articles [8,9] showed that entrapped or escaping air cause cracking near the surface which results in a damage of the structure.Cracks appear between the entrapped air areas and outer boundaries.While such phenomena were not investigated in this contribution it is believed that an explicit taking into account of the air flow in the simulation is a first step to their modelling.

Fig. 4 .
Fig. 4. Example 1: distribution of seepage flux at the inner slope of the dike.

Fig. 6 .
Fig. 6.Example 2: evolution of total volume of water in dike.