Air trapping problem during infiltration on the large areas

The process of flow modeling in unsaturated porous medium is often found in many fields of sciences: geology, fluid mechanics, thermodynamics, microbiology or chemistry. Problem is relatively complicated due to complexity of the system which contains three phases: water, air and soil skeleton. The flow of water in such a medium can be described using two-phase (2PH) flow formulation, which accounts the inflow of air and water phases, or with simplified model known as Richards (RE) equation where only water flow is taken into account. In many well known programs available in the market (like SeepW, STOMP) the primary interest is only the water flow and the flow of air is omitted. As a result Richard equation in used more often. It’s main assumption is that pore air is continuous and has connection with atmospheric air which is equivalent to infinite mobility of the air phase during all simulation. 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 infiltration process occurring in the large areas. An irrigation effect of rice fields with shallow water table has been investigated. To assess the impact of the gas phase various lengths of the infiltration zone have been considered. 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.


Introduction
Infiltration of water in saturated porous medium has been studied especially through the Darcy equation (single phase flow). For unsaturated soils this formula was extended to the Richards equation. Many programs available in the market are designed to compute problems only in the case of single phase flow (usually water) with assumption that gas phase does not affect the infiltration of liquid phase into soil. This is caused mainly by mobility of gas phase which is infinite -principle assumption of Richard equation. However mobility of air phase is not infinite but 50-60 times bigger than mobility of water phase. Experiment performed by [1], showed that agreement between Richard equation and two phase flow model can be obtained only if we assume that viscosity of air is about 100 times bigger than viscosity of water. So, even if gas phase present in soil has a free connection with atmospheric air we will obtain discrepancies between both models. Another important feature is relative conductivity. Additional decrease in air mobility we can obtain when conductivity of water is definitely bigger than for air. This problem occurs when saturation is very close to its final value [2].
Mobility of air phase can be also limited by various boundary conditions: heavy rainfall [3,4], surge flooding [5,6], highly impermeable inclusions [7,8] and irrigation of agricultural areas [9]. Problem is especially important when water table level is close to the ground surface and the infiltration takes place on large area. This causes an increase in pore air pressure which decreases the infiltration rate. In addition, huge amounts of air are partially trapped between water table level and the water that infiltrates from the surface of the soil. The consequences of this situation have been evaluated theoretically and in laboratory experiments for example [10]. The most important effect is closing and compressing the air, which significantly slows down the infiltration [11]. Another possibility is unstable flow and fingering process caused by air bubbles that escape through the domain to the atmosphere. Otherwise just a little is known about the phenomenon of air trapping under crops [9]. Effect is of particular interest in flooded rice cropping in arid regions. This article is a continuation of article [9] where one of the conclusions was comparison of two solutions -Richards equation and two-phase flow model.

Equations
Single phase flow occurs when the pores are completely filled with one fluid phase, which is the context of this article is water. Therefore saturation and relative water permeability are equal one. Under these assumptions governing equation can be written in the following form: where  is phase occupying the pore space,  -density, nporosity, t -time,  -intrinsic permeability, μ -dynamic viscosity, ppressure in fluid, g -gravity, z -depth.
Unfortunately equation (1) can describe only water flow in saturated soils. While for unsaturated is necessary to add saturation and relative permeability components: this form of equation is known as Richards equation. It is used to solve problems related to the flow of water in many popular programs like STOMP or SeepW. If we want to consider also the role of the gas phase it is necessary to solve the system of two equations -one for water and one for air. Both the two phase flow model (2PH) and Richards equation (RE) are difficult to solve numerically because of nonlinearities describing the relations between e.g. saturation and permeability or capillary pressure and type of the soil (clay, sand, silt).
Mathematical model used in this article is based on the equation (2) and 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. 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 and pressure difference between the non-wetting fluid (air) and the wetting fluid (water) is called capillary pressure or suction: 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r is residual saturation of fluid, Sew -effective water saturation, pg -scaling pressure in the retention function, ng and mg are the exponents in the retention function. Permeability of the medium describing the Mualem-van Genuchten model: Besides definition of the constitutive relationship it is also necessary to select appropriate par of unknowns. Three most popular are: pressure formulation and pressuresaturation formulation. Unfortunately pressure formulation has the big disadvantage that, for the solution of the system of equations, the capillary pressure gradient have to be greater than 0. In the case of infiltration or heterogeneities this gradient is very small or even 0. So, this restriction makes the use of pressure formulation impossible. Therefore we use the mixed formulation (pressure-saturation) due to possibility of disappearing one of the phases (such as the air phase during variably saturated flow in soils) [12,13].

Numerical implementation
Numerical solutions for Richards equation and 2PH model were obtained using an in-hose code developed by A. Szymkiewicz. It is based on the finite volume method (FVM) for the spatial discretization combined with the first-order fully implicit temporal discretization scheme. Detailed description of the code can be found in [12]. It is important to note that the same scheme was used for both solutions in order to reduce the influence of the numerical errors. Spatial discretization is done using the vertex-centred finite volume method. In this approach the spatial domain is firstly covered with a finite element grid and after that a dual grid is created by defining a finite volume around each node.
In the case of Richards model for each node we will obtain one equation which describe the water flow. While for 2PH model we will get two -for water and air respectively. Thus, a system of nonlinear algebraic equations is obtained, which allow computing the values of the primary unknowns at the new time level based of the knowledge of the values from the old time level. Due to the stability of the solution as the pair of primary unknowns have been chosen water pressure and water saturation. Other values such as relative permeability, capillary pressure or air pressure are calculated on the basis of this to pair of unknowns. The nonlinear algebraic system arising at each time step is solved using Newton-Raphson iterative scheme.

Geometry, material parameters and boundary conditions
The influence of air phase on the water flow during infiltration in large area is investigated for three different length of the irrigation zone. The first test represents situation in which water is infiltrating with 10 m. Two further simulations show what will happen when irrigation length is increased or decreased by 5 ( fig. 1). We consider homogeneous area placed on an impermeable base. Water table level is placed 2 m bellow the ground surface. Soil parameters correspond to loamy sand and were taken from [14] and are as follows: Swr = 0.16, Sar = 0.03, pg = 1308 Pa, ng = 1.89, mg = 0.471. For both fluid phases it is necessary to implement density and viscosity which are respectively: w=1000 kg/m 3 , μw=10 -3 Pa·s for water and a=1.22 kg/m 3 , μw=1.57·10 -5 Pa·s for air.

Fig. 1. Dike geometry and boundary conditions used in simulations.
Initially area above the water table is in unsaturated state. Water pressure in the domain is assumed to varying hydrostatically form 9.81 kPa on the bottom to -19.62 kPa on the top. Irrigation stream has been modelled as water table keeping level of 0.5 m. Air pressure in the whole area was assumed to be constant and equal to 0 (atmospheric value). Left site of the domain is impermeable for water and air. While on the right site flow of both phases is possible. Remain part of the upper boundary where there is no irrigation is considered as a permeable for air.
Numerical grid was performed using quadrilateral elements with computational points in the grid nodes (vertex centred method). Here we show results for a 6000 elements and 6231 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. 2-3 and they confirm suppositions put forward by [Hammecker, 2003]. Figure 2 shows the distribution of water saturation in the areas at a mid stage and near the end of the simulation. It can be see that when we will use Richards equation an area between water tables is fully saturated at 4800 s. While in 2PH model there is a large amount of air trapped between infiltrating stream and ground water table. Over time, the air volume in the area decreases. Some part of the gas phase escape the domain in the form of bubbles, rest is dissolved in water. As a result saturation maps of both models are more similar to each other. Despite this, in the end of the simulation, mass of water present in the domain is different for both models (fig. 3). This proves that some part of air did not dissolve and still exit in the area.  Fig. 3 shows the distribution of the water mass for three different lengths of the irrigation stream. With the increase of the infiltration zone differences between Richards equation and 2PH model are also increase. It is caused by isolating bigger and bigger air areas and the gas phase compression, which is taking into account only by the 2PH model. While in Richards, water phase can infiltrate freely without necessity to displace the air phase. This is related to the main principle of the Richards model which assumes that air present in the soil has an infinite mobility and can connect freely with atmospheric air. To better see the inflow of mobility on the obtained results a few additional simulations has been done. Fig. 4 presents a water saturation values for different values of the viscosity which are varying from 60 to 1000. It is clear to see that with increasing values of the viscosity the solution obtained from 2PH model is closer to RE but agreement between these two solutions is achieved for viscosity of air about 1000 times smaller than viscosity of water. This result is not in agreement with [1] where author noticed that 2PH model is giving the same results as RE for viscosities differing by 100.
For each simulation, air and water pressure in point A ( fig. 5a and 4b) were also checked. In Richards equation saturation is the same for model with 10 and 15 m infiltrating stream. While using the 2PH model we can see that achieving the appropriate water pressure occurs definitely later and depends also on the length of the irrigation zone. For 5 m infiltration both pressures are less compared to the other values gained for 10 and 15 m. The reason is relatively small irrigation zone which connects with water table near the point A. In the other cases this place is far away from the analyzed point so the pressure rise process is lengthened.

Conclusions
Air trapping problem presented in this article is a continuation of [9]. Large agricultural areas with shallow water table are very interesting in the case of air trapping especially during irrigation process. Numerical simulation presented in the article showed huge differences between the unsaturated flow model (Richards equation) and the two-phase model which is in agreement with conclusions described by [9]. Differences between models increase parallel with infiltrating zone. It is related to the area from which the air present in the soil can freely escape to the atmosphere. In Richards equation is not possible to include this effect due to assuming the infinite mobility of air phase. While in 2PH model air has a finite mobility and infiltrating water needs to displace the gas phase during the whole process. In work presented by [1] mentioned that differences between models resulting from a finite mobility of air can be solved when viscosity of air will be about 100 times bigger than viscosity of water. However, curves on fig. 4 showed that this is not true even if viscosity difference is much bigger than 100. Similar results were obtained only for values differing by 1000. This can be caused by relatively big domain and boundary condition which hampering the connection between air present in the soil and atmospheric air.