Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2011, Article ID 163020, 23 pages
Research Article

Using Global Characteristics of a Centrifuge Outflow Experiment to Determine Unsaturated Soil Parameters

1Faculty of Mathematics, Physics and Informatics, Comenius University Bratislava, Mlynská dolina, 84248 Bratislava, Slovakia
2Department of Mathematical Analysis, Research Group NaM2, Ghent University, Galglaan 2, B 9000 Gent, Belgium

Received 29 March 2011; Revised 8 September 2011; Accepted 15 September 2011

Academic Editor: J. Rodellar

Copyright © 2011 Jozef Kačur et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Several centrifugation scenarios enabling the determination of soil parameters for saturated-unsaturated flow in porous media are presented, investigated, and discussed. Only global characteristics of the infiltration process in a sample are used, so that only simple, noninvasive measurements are performed. The characteristics can be transient measurements of the rotational momentum, or of the gravitational center, or of the water amount injected and expelled from the sample. No information about the saturation or the head distribution in the sample is required. This setup is different from the common multioutflow experiments. We give numerical proof that this method allows for fast determination of soil parameters in comparison to traditional measurements based on equilibrium conditions. The mathematical model of infiltration is represented by Richards' strongly nonlinear and degenerate equation expressed in terms of soil parameters in the van Genuchten-Mualem ansatz. The parameter identification process is realized in an iterative way applying the Levenberg-Marquardt method. Numerical experiments support the efficiency of the analyzed method and allow one to identify the optimal centrifugation scenario for imbibition and drainage to be applied when using global characteristics.

1. Introduction

To predict flow and solute transport in soils concerns the soil hydraulic properties in terms of soil parameters, which are required in the governing mathematical model. This model is expressed in terms of saturation and pressure head in Richards’ equation, which is a nonlinear and degenerate parabolic equation with free boundaries between saturated and partially saturated zones and between dry and partially saturated zones.

The capillary pressure and hydraulic permeability functions linking the saturation and pressure head are expressed using the van Genuchten-Mualem ansatz by means of soil parameters. To measure these soil parameters is usually time consuming and expensive, specifically for low-conductive materials. The determination can be accelerated by using centrifugation. The measurement of capillary pressure curves with a centrifuge was initiated in [1]. With the development of more efficient numerical tools, this method became more popular in the last decades. A more detailed overview on this topic can be found in [2] (see also citations there).

Recently, the method of centrifugation has been applied in [3, 4] where the equilibrium and transient analysis at a set of rotational speeds has been used for the determination of soil parameters. In [4], the distribution of the saturation in equilibria (linked with the corresponding rotational speeds) is measured via electrical signals from electrodes installed in the sample.

The majority of applications of centrifugation are based on the measurements in a series of equilibria reached in the sample by application of a fixed rotational speed and by control of the boundary conditions. This can require a very long time of centrifugation and requires information on the error, since reaching equilibrium is an asymptotic process. A first approach applies the series of equilibria and the amount of expelled water to obtain the retention curve, see overview in [5]. This approach is slow however. A second approach requires transient data from inside the sample at one or more points, which results in relatively expensive measurements. This approach is similar to the common multistep outflow experiments to determine the soil retention curve. For example, [6] shows that such an approach combined with an inverse method allows one to determine hysteretic hydraulic properties of the soil retention curve.

Our main goal is to replace the internal measurements with global characteristics such as rotational momentum of the water infiltrated into the sample, gravitational center of the water content, and time evolution of the amount of water going into and coming out of the sample. It is possible to measure these variables during rotation of the centrifuge. Hence, the considered measurements do not include pointwise saturation or head distributions over the sample. Doing so, however, adds some constraints to the approximation method used: conservation of mass and correct prediction of the interfaces (wetness and saturated fronts) will be important. This is because they have a large influence on the computed rotational momentum and gravitational center, meaning that errors would have a large influence. We do not focus here on how these global measurements can be done practically. Based on the results we present here, a project has started with that specific goal, testing various measuring techniques in a custom-built, table-top centrifuge.

The mathematical model based on Richards’ equation is highly nonlinear with free boundaries and a sharp front on the interface between partially saturated and dry zones in the sample. On the other hand, the free boundary between the fully saturated and the partially saturated front is difficult to identify (when they are expressed in terms of saturation, see Figure 2). It is therefore a difficult task for numerical approximations to supply experimenters with accurate values to use in the determination of model parameters. Additionally, the computation has to be able to simulate long centrifugation runs (sometimes weeks). The second goal in this paper is therefore to present a new approximation method which is accurate and efficient and which can supply us with the needed global characteristics.

Using only global characteristics makes the inverse determination process of the soil parameters more difficult, so we investigate the sensitivity of the global characteristics, which could be small. Numerical experiments previously done by us indicate that a series of equilibria are not suitable to obtain strong sensitivity of global characteristics on the soil parameters. Moreover, in this setup, the centrifugation has to continue for a very long time to obtain sufficient experimental data. As an extra complication, creating equilibria by centrifugation (with zero input and output flux) means that some regions in the sample exhibit drainage while others undergo imbibition. Therefore some aspects of hysteresis should then be taken into consideration, but this increases the number of unknown parameters considerably. The dynamics behind the hysteresis phenomenon of the retention curve, and the optimal way of modeling it, is still being discussed in the literature. Practical numerical modeling has been realized only under strict conditions, see, for example, [68]. For this reason, the centrifugation scenarios we use are designed in such a way as to avoid the hysteresis phenomenon, using only the extremal drainage curve and imbibition curve. In a previous article [9], only drainage is considered in a start-stop scenario. Here we extend this method and present a scenario for drainage and imbibition that allows continuous operation of the centrifuge.

The centrifugation scenario for imbibition consists of a water chamber on top of the sample that allows infiltration into an originally dry sample. On the other hand, the soil parameters in drainage mode can be obtained by centrifugating a fully saturated sample with zero input flux. In both cases it is assumed that no compactification occurs. Injecting water (by centrifugation) into an originally dry sample will give rise to a fully saturated zone in the front of the sample followed by a partially saturated and a dry zone. Tracking the water movement in this system leads to global characteristics that are sensitive to the soil parameters. As determination of the saturated hydraulic permeability (𝐾𝑠) is well known, see, for example [2], we can reduce the number of soil parameters by assuming 𝐾𝑠 to be known.

The presented mathematical model and its numerical realization can be applied in many different centrifugation scenarios. The accurate numerical approximation is based on the following principles:(i)application of a mathematical model for the evolution of the wetness front (interface between partially saturated and dry region), which was developed in [3], see Section 2;(ii)application of moving grid points which significantly increase the numerical accuracy and effectiveness;(iii)numerical modeling of the free boundary between saturated and partially saturated zones, based on global water mass balance, extending the method developed in [9];(iv)expressing the governing Richards’ equation or in terms of saturation, or in terms of pressure head, depending on the distance to both free boundaries which separate the partially saturated zone from the dry and fully saturated zones;(v)approximation of the governing mathematical model (1D) by a system of ODE, the so-called method of lines (MOL).

In Section 2 the mathematical model is introduced, and in Section 3 the numerical method is used. Numerical experiments are discussed in Section 4. In Experiment 1 we present the solution of infiltration into the sample; Experiment 2 discusses the applicability of the method for the determination of soil parameters in imbibition mode, Experiment 3 presents the determination of the saturated conductivity, and Experiment 4 contains the determination of the soil parameters in drainage mode. Finally, in Experiment 5 we demonstrate the applicability of the method for a larger scale of soil parameters, as in the first experiments only one set of “standard soil parameters’’ (corresponding to layered clay) has been used. Moreover, we discuss the infiltration of water into the sample during the preparation period. The determination of soil parameters is realized iteratively using the Levenberg-Marquardt method (LM), as in [6]. In all these experiments we demonstrate the practical applicability of global characteristics.

2. Mathematical Model

We consider a one-dimensional dry sample in the form of a tube which starts (left boundary) at the distance 𝑟=𝑟0 from the center of the centrifuge and ends at the distance 𝑟=𝑟0+𝐿, see Figure 1. A water chamber is placed at position 𝑟(𝑟00,𝑟0)of the centrifuge. When the dry sample, defined as a sample from which no longer water can be removed, comes in contact with water, different zones will arise. First, a partially saturated zone will appear, characterized by capillary flow. In this zone, air and water are present in the pores. Secondly, a saturated zone will appear. This zone is characterized by pores which are completely filled by water.

Figure 1: Schematic representation of the centrifuge. A water chamber to the left of height , sample for 𝑟(𝑟0,𝑟0+𝐿), and outflow chamber to the right. Full saturation up to 𝑟0+𝑠1 and dry zone to the right of 𝑟0+𝑠2.
Figure 2: Time evolution of the saturation, 𝑡=𝑖𝑇out/20,𝑖=0,,20;𝑇out=5016. Left curve at 𝑡=0, right curve for 𝑡=𝑇out.

Saturated flow (piezometric head 0) in porous media under centrifugation is modeled by Darcy’s equation, while the unsaturated part (<0) is governed by Richards’ equation as follows:𝐾𝑠𝜕𝑟𝜕𝑟𝜔2𝑔𝑟𝜕=0,0,𝑡𝜃=𝐾𝑠𝜕𝑟𝜕𝑘(𝜃)𝑟𝜔2𝑔𝑟,<0,(2.1) where 𝜃 is the saturation of the porous media, 𝜔 the angular speed of rotation (in radians per second), 𝐾𝑠 the saturated hydraulic conductivity, 𝑔 the gravitational constant, and function 𝐾𝑠𝑘(𝜃) the hydraulic conductivity in the unsaturated region. Denote by 𝑢=(𝜃𝜃𝑟)/(𝜃𝑠𝜃𝑟) the effective saturation, where 𝜃𝑠 is the volumetric water content at saturation and 𝜃𝑟 the residual volumetric water content. We have 𝑢(0,1), since 𝜃(𝜃𝑟,𝜃𝑠). We consider soil hydraulic properties which were proposed in [10]1𝑢=1+(𝛾)𝑛𝑚,(,0),𝑘(𝑢)=𝑢1/211𝑢1/𝑚𝑚2,(2.2) where 𝑚=11/𝑛, 𝑛>1, and 𝛾=(21/𝑚1)1𝑚/𝑏 are empirical soil parameters and 𝑏 is the bubbling pressure. Hence, flow in an unsaturated region can be rewritten in terms of saturation as𝜕𝑡𝑢=𝐾𝑠𝜕𝑟𝐷(𝑢)𝜕𝑟𝜔𝑢2𝑔𝑘(𝑢)𝑟,(2.3) where1𝐷(𝑢)=𝜃(𝑛1)𝛾𝑠𝜃𝑟×𝑢1/21/𝑚1𝑢1/𝑚𝑚11𝑢1/𝑚𝑚2.(2.4) Equation (2.3) is strongly nonlinear and degenerate (𝐷(0)=0, 𝐷(1)=). As a consequence of this, the propagation of the wetness front into the dry region will proceed with finite speed (i.e., there appears an interface or free boundary). The first equation in (2.1) can be integrated, which leads to𝜔(𝑟)=2𝑟2𝑔2+𝐶1𝑟+𝐶2.(2.5) During the centrifugation, the water from the reservoir is pushed into the dry region of the sample and creates a saturated and unsaturated subregion in the sample. Denote by 𝑠1(𝑡)(0,𝐿) the free boundary which separates the saturated and unsaturated regions in the sample, and denote by 𝑠2(𝑡) the free boundary separating the partially saturated zone from the dry zone (wetness front), 𝑠1(𝑡)<𝑠2(𝑡)<𝐿. To obtain the solution for 𝑟(𝑟0,𝑟0+𝑠1(𝑡)), the 2 integration parameters in (2.5) need to be determined from boundary conditions. We can easily determine that𝑟0=𝜔,𝑡22𝑔(𝑡)2𝑟0𝑠(𝑡),1(𝑡),𝑡=0,(2.6) since the piezometric head at 𝑟0 is equal to the pressure of the water in the reservoir (contained in the interval (𝑟0(𝑡),𝑟0)) and the piezometric head on the interface 𝑟0+𝑠1(𝑡) is zero. After elimination of the integration parameters we obtain the solution of (2.1) in the saturated region 𝑟(𝑟0,𝑟0+𝑠1(𝑡)) in the form𝜔(𝑟,𝑡)=22𝑔𝑟𝑟02𝑟𝑟01𝑠12𝑟0+𝑠1+2𝑟0.(2.7) From this we can derive the flux 𝑞(𝑡) at the interface as𝑞(𝑡)=𝐾𝑠𝜕𝑦𝑟0||+𝑦𝑦=𝑠1(𝑡)𝜔2𝑔𝑟0+𝑠1(𝑡)=𝐾𝑠𝜔212𝑔𝑠1(𝑡)2𝑟0𝑠1(𝑡)+𝑠1(𝑡)2+(𝑡)2𝑟0,(𝑡)(2.8) where (𝑡) and 𝑠1(𝑡) have to be determined. Since the flux of water along (𝑟0,𝑟0+𝑠1(𝑡)) is constant (at fixed 𝑡) due to mass balance, it follows thaṫ(𝑡)=𝑞(𝑡)=𝑓1,𝑠1,(2.9) where 𝑞 is obtained from (2.8). We proceed by constructing a model for ̇𝑠1(𝑡). Because of mass balance reasons, the flux 𝑞(𝑡) is used to increase the saturated and unsaturated subregions. The crucial point is to obtain the flux 𝑞int(𝑡) entering the unsaturated subregion. The unsaturated region is placed between two free boundaries 𝑟(𝑟0+𝑠1(𝑡),𝑟0+𝑠2(𝑡)). We note that 𝑢(𝑠2(𝑡),𝑡)=0 and 𝑢(𝑠1(𝑡),𝑡)=1, forall𝑡. Then 𝑢 in the unsaturated region is governed by (2.1) with the moving domain of solution (𝑠1(𝑡),𝑠2(𝑡)) (we shift the original domain by 𝑟0). We transform the moving domain to a fixed domain, using the transformation 𝑦=(𝑟𝑠1(𝑡))/(𝑠2(𝑡)𝑠1(𝑡)). Denote 𝑢(𝑦,𝑡)=𝑢(𝑟,𝑡), and, for simplicity, in the following we drop the overline, writing shorthand 𝑢 for 𝑢(𝑦,𝑡). We have that𝜕𝑡𝐾𝑢(𝑦,𝑡)=𝑠𝑠21(𝑡)2𝜕𝑦𝐷(𝑢)𝜕𝑦𝑢𝑠21𝜔(𝑡)𝑘(𝑢)2𝑔𝑟0+𝑠1(𝑡)+𝑦𝑠21+(𝑡)̇𝑠1(𝑡)(1𝑦)+̇𝑠21(𝑡)𝑦𝑠21𝜕(𝑡)𝑦𝑢,(2.10) where 𝑠21(𝑡)=𝑠2(𝑡)𝑠1(𝑡), and with boundary and initial conditions𝑢(0,𝑡)=1,𝑢(1,𝑡)=0;𝑢(𝑦,0)=𝑢0(𝑦).(2.11) Modeling the interface 𝑠2(𝑡), we follow [3] where problems similar to (2.10) have been studied. We compute first the order of degeneracy of 𝐷:lim𝑧0𝐷(𝑧)𝑧𝑝=1𝜃(𝑛1)𝛾𝑠𝜃𝑟𝑚2,(2.12) where 𝑝=1/2+1/𝑚. Therefore,̇𝑠2𝐾(𝑡)=𝑠𝑚2𝑝𝜃𝑠𝜃𝑟1𝛾(𝑛1)𝑠21𝜕(𝑡)𝑦𝑢(1,𝑡)𝑝1,where𝑝=2+1𝑚.(2.13)

Now, the solution of (2.10)–(2.13) defines a flux 𝑞int(𝑡) which is the flux of the solution at 𝑦=0. Finally, we close our system with the mass balance condition 𝑞(𝑡)=̇𝑠1(𝑡)+𝑞int(𝑡), from which we obtaiṅ𝑠1(𝑡)=𝑞(𝑡)𝑞int(𝑡).(2.14) The final system (2.9)–(2.14) can be solved by approximating it with a corresponding ODE system.

This model has to be modified when the wetting front (𝑠2) reaches the right boundary. Note that we do not consider the case where the injection becomes empty (=0), as in this scenario we want to avoid drainage (which starts from the left when =0). If however at some 𝑡=𝑇out we have >0 and 𝑠2=𝐿, the dry zone is no longer present. At that instance, we start to control the amount of water expelled from the right boundary of the sample. Experimentally, the water level 𝑤out in the water-collecting chamber situated on the right of the sample can now be measured. In the mathematical model we obtain it from the flux on the boundary which is linked with ̇𝑤out and with the global water mass balance.

In the case of a drainage scenario, the initial condition is a saturated sample with no injection water chamber. Only drainage will occur, and the outflow can be controlled in the same way as in the above scenario when the wetness front reaches the outflow boundary. The difference with [9] lies in the attempt to have a continuous operation of the centrifuge. In [9] it was shown that the soil parameters can be recovered from global characteristics when starting from a partially saturated sample and applying first a period with free outflow, then a period with no outflow up to equilibrium, and finally a period with again free outflow. The drawback of that approach is firstly that the operation of the centrifuge is not continuous and secondly that during the no-outflow period imbibition takes place in part of the sample. In this paper we show in the experiments that starting from a fully saturated sample, the soil parameters can also be retrieved, with some caveats however.

Remark 2.1. Practically, the ground sample will have to be contained in the centrifuge via two fixed filters (porous stone) at the start and finish. For simplicity this is neglected in this study, but a realistic model will have to take them into account by using variable groundwater characteristics like porosity and saturated conductivity, ….

3. Numerical Method

The 1D mathematical model results in a coupled system of PDE and ODE (2.9)–(2.14) and is therefore a good candidate for applying first a space discretization and then solving the resulting ODE system (MOL method). The space discretization can reflect the expected profile of the solution in 𝑥(𝑠1(𝑡),𝑠2(𝑡)) which was transformed to a fixed domain 𝑦(0,1). Apply the space discretization0=𝑦0<𝑦1<<𝑦𝑖<𝑦𝑁=1,(3.1) and define 𝛼0=0,𝛼𝑖=𝑦𝑖𝑦𝑖1,𝑖=1,,𝑁. This discretization corresponds to moving grid points 𝑥𝑖(𝑡)=𝑟0+𝑠1(𝑡)+𝑦𝑖𝑠21(𝑡) in the sample. Since the solution has a sharp front at the wetness interface (see Figure 3), an increasing density of grid points near the boundary point 𝑦=1 is chosen. The used transformation then guarantees a good approximation accuracy of space derivatives at the sharp front during the entire computation. For the space discretization equation (2.10) is integrated over the domain 𝐼𝑖=(𝑦𝑖𝛼𝑖/2,𝑦𝑖+𝛼𝑖+1/2),forall𝑖=1,,𝑁1. Denote 𝑢𝑖(𝑡)𝑢(𝑦𝑖,𝑡) and 𝜕𝑡𝑢(𝑦,𝑡)̇𝑢𝑖(𝑡) in the interval 𝐼𝑖, forall𝑖=1,,𝑁1. Approximate𝜕𝑦𝑢||𝑦=𝑦𝑖+1/2𝑢𝑖+1(𝑡)𝑢𝑖(𝑡)𝛼𝑖+1=𝜕+𝑢𝑖,(3.2) where 𝑦𝑖+1/2=𝑦𝑖+𝛼𝑖+1/2, and similarly approximate 𝜕𝑦𝑢|𝑦=𝑦𝑖1/2 and denote it by 𝜕𝑢𝑖. Also denote 𝑢𝑖+1/2=(𝑢𝑖+1+𝑢𝑖)/2 and 𝑘𝑖+1/2=𝑘(𝑢𝑖+1/2) (similarly for 𝑖1/2). Then, the approximation of (2.10) at the point 𝑦=𝑦𝑖, 𝑖=1,,𝑁1, reads as follows:̇𝑢𝑖=2𝐾𝑠𝛼𝑖+𝛼𝑖+1𝜃𝑠𝜃𝑟1𝑠221×𝐷𝑖+1/2𝜕+𝑢𝑖𝐷𝑖1/2𝜕𝑢𝑖𝜔2𝑠21𝑔×𝑘𝑖+1/2𝑟0+𝑠1+𝑦𝑖+1/2𝑠21𝑘𝑖1/2𝑟0+𝑠1+𝑦𝑖1/2𝑠21+̇𝑠11𝑦𝑖+̇𝑠2𝑦𝑖𝑑𝑧;𝑦𝑖||||𝑑𝑧𝑧=𝑦𝑖,̇(3.3)(𝑡)=𝑓1,𝑠1,(3.4)̇𝑠1(𝑡)=𝑓1,𝑠1𝑓2𝑠1,𝑠2,𝑢1,𝑢2,(3.5)̇𝑠2(𝑡)=𝑓3𝑠1,𝑠2,𝑢𝑁2,𝑢𝑁1,(3.6) where (𝑧;𝑦𝑖) is the Lagrange polynomial of the second order crossing the points (𝑦𝑖1,𝑢𝑖1(𝑡)), (𝑦𝑖,𝑢𝑖(𝑡)), and(𝑦𝑖+1,𝑢𝑖+1(𝑡)), and function 𝑓3 is given by𝑓3𝑠1,𝑠2,𝑢𝑁2,𝑢𝑁1𝐾=𝑠𝑚2𝑝𝜃𝑠𝜃𝑟1𝛾(𝑛1)𝑠21(𝑡)𝑑𝑧;𝑦𝑁||||𝑑𝑧𝑧=1𝑝,(3.7) where 𝑝=1/2+1/𝑚 and (𝑧;𝑦𝑁) is the Lagrange polynomial crossing the points (𝑦𝑁2,𝑢𝑁2(𝑡)), (𝑦𝑁1,𝑢𝑁1(𝑡)), and (1,0). The function 𝑓2 represents the numerical approximation of the flux 𝑞int and reads as follows:𝑓2𝑠1,𝑠2,𝑢1,𝑢2=𝐾𝑠𝑠21𝜃𝑠𝜃𝑟𝑑+𝑧;𝑦0||||𝑑𝑧𝑧=0𝜔2𝑔𝑠1,(3.8) where +(𝑧;𝑦0) is the Lagrange polynomial crossing the points (0,0), (𝑦1,1), and (𝑦2,2), with𝑖1=𝛾1+𝑢𝑖1/𝑚1/𝑛,𝑖=1,2.(3.9) We note that the flux at 𝑦=0 must be expressed in terms of head (instead of saturation) since at 𝑦=0 we have 𝑢0=1, 0=0, and 𝐷(1)=+. We replace ̇𝑠1, ̇𝑠2 in (3.3) by 𝑓1𝑓2 and 𝑓3, respectively, to obtain the ODE systeṁ𝑢𝑤=𝐹(𝑡,𝑤),𝑤=1,,𝑢𝑁1,,𝑠1,𝑠2,(3.10) and solve it by an ODE solver for stiff systems. This system must be completed by the initial state 𝑤(0)=[𝑢1(0),,𝑢𝑁1(0),(0),𝑠1(0),𝑠2(0)]. Due to the first experimental scenario (infiltration into a dry region), all components of 𝑢 would be zero, and also 𝑠1(0)=𝑠2(0)=0 and (0)=0. This would be a singular case for our ODE system (initial and boundary conditions are not compatible), and hence, a slight regularization is applied, in the sense that we take 𝑠2(0)>𝑠1(0)>0, with 0<𝑢<1 in between.

Figure 3: Time evolution of the head, 𝑡=𝑖𝑇out/20,𝑖=0,,20;𝑇out=5016. Left curve at 𝑡=0, right curve for 𝑡=𝑇out.

A problem in the above approach to discretizing Richards’ equation is that typically mass balance is not well kept, see, for example, [11], which is problematic as errors in mass balance will have a big influence on the global characteristics. Specifically to this approach is also that 𝑠1 is not accurately determined. We do not follow the solution of [11] which consists of using a mass-pressure scheme. Instead, like in [9], we drop (3.5) and replace it by the global mass balance equation𝜃𝑠𝜃𝑟𝑠1(𝑡)+𝑁𝑖=0𝑐𝑖𝑢𝑖+𝑤out(𝑡)=(𝑡)(0),(3.11) where𝑐𝑖=𝑠21𝛼(𝑡)𝑖+1+𝛼𝑖2𝑐,𝑖=1,,𝑁1,0=𝑠21𝛼(𝑡)02,𝑐𝑁=𝑠21𝛼(𝑡)02.(3.12) This mathematical model can be used up to the time 𝑡=𝑇out when 𝑠2(𝑇out)=𝐿. At that moment it is needed to replace 𝑠2 by 𝑢𝑁𝑢(1,𝑡) and put 𝑠2𝐿; that is, we drop (3.6) and add the ODE for ̇𝑢𝑁 to (3.3) (created on the same argumentation as for 1𝑖𝑁1). Moreover, the ODE system is completed bẏ𝑤out(𝑡)=𝐾𝑠𝐷𝑢𝑁𝜕𝑦𝑟0||+𝑦𝑦=1𝜔2𝑔𝑘𝑢𝑁𝑟0+𝑠1(𝑡)(3.13) linking the outflow flux and 𝑢𝑁. Now, the approximation consists of 𝑁+2 ODE and the algebraic equation (3.11). This system can be written in the form𝑀(𝑡,𝑧)̇𝑧(𝑡)=𝑓(𝑡,𝑧),(3.14) where𝑧=1,,𝑖0,𝑢𝑖0+1,,𝑢𝑁1,𝑠1,𝑢𝑁,,𝑤out.(3.15) In this, the index 𝑖0 is chosen suitably large so that towards the saturated part pressure head is used as unknown, while towards the dry part the saturation is used. Practically, for the computations, we have chosen 𝑖0=𝑁/2. Note also that the 𝑁th row of the matrix 𝑀 consists of zeros as it corresponds to the algebraic equation modeling mass balance.

3.1. Centrifugation of an Initially Saturated Sample with Outflow Control

We shortly discuss the second scenario where the sample undergoes drainage during the centrifugation. Now, the input flux is zero and there is an ODE for 𝑢0. The governing ODE system is obtained from (3.3) with 𝑖=0,,𝑁1, where 𝑠1=0, 𝑠2=𝐿, and 𝑠21=𝐿, =0. An ODE is added for ̇𝑢0 similarly as for ̇𝑢𝑁 before. The outflow 𝑤out is governed by (3.13), where 𝑟0+𝑠1(𝑡) is replaced with 𝑟0+𝐿. The system is completed with the mass balance equation (3.11), where 𝑠1=0 and (𝑡)(0) is replaced by 𝑀𝑤 which represents the mass of water in the fully saturated sample. The system (3.14) has to be rewritten in terms of𝑧=0,,𝑖0,𝑢𝑖0+1,,𝑢𝑁,,𝑤out.(3.16)

3.2. Centrifugation of a Saturated Sample with Input and Output Control

In order to allow one single experimental run to determine the soil parameters including the saturated conductivity 𝐾𝑠, we discuss how centrifugation of a fully saturated sample with injection chamber and outflow fits into the presented modeling. One can start with a saturated sample and a filled injection chamber and measure 𝐾𝑠. The right boundary of the sample is free. In this case, the injected water amount is equal to the amount of output water. This significantly changes the measured rotational moment or gravitational center and consequently increases the reliability of the determination of 𝐾𝑠. In this case there is one governing ODE in terms of , given bẏ(𝑡)=𝐾𝑠𝜔2𝐿2𝑔𝐿2(𝑡)2+2𝑟0,(𝐿+(𝑡))(3.17) and we can proceed centrifugation up to =0.

4. Numerical Experiments and Discussion

As only numerical experiments are performed, variables that have correct order of magnitude are chosen, dropping the units. Unless clearly noted otherwise, the following values are used: 𝑟0=30, 𝐿=10, 𝜔=20, 𝐾𝑠=2.4105, 𝜃𝑟=0.02, 𝜃𝑠=0.4, 𝛾=0.0189, and 𝑛=2.81. For the space discretization 𝑁=40 grid points are used with geometrical distribution: a first space interval with width 𝑑1=𝐿/20, and then 𝑑𝑖+1=𝑞𝑑𝑖 (𝑖=1,,39) with 𝑞<1. Furthermore, the surface area of the sample and the water chambers is assumed to be 1. The chamber for expelled water is situated in 𝑟(𝐿,𝐿+𝑅), with 𝑅=4. At the start it is assumed that no water is present in the outflow chamber (which remains the case up to 𝑡<𝑇out). In the case of a filled injection chamber, initial water level 0=6 is taken.

The global characteristics that we assume can be measured accurately are the rotational momentum 𝑀𝑟 (rotational kinetic energy of the water mass, both in the sample and chambers), gravitational center 𝐺, and water mass in the injection and outflow chamber, respectively, and 𝑤out (as the surface is 1). The partially saturated zone was transformed to the domain 𝑦(0,1). The contribution of the mobile water in the partially saturated zone to the global characteristics is denoted by superindex (𝑝) and can be calculated from𝑀𝑟(𝑝)𝑠(𝑡)=21(𝑡)𝜔2210𝑟0+𝑠1+𝑠21(𝑡)𝑦2𝑀𝑢(𝑦,𝑡)d𝑦,𝑤(𝑝)(𝑡)=𝑠21(𝑡)10𝐺𝑢(𝑦,𝑡)d𝑦,(𝑝)(𝑡)=𝑠21(𝑡)10𝑟0+𝑠1+𝑠21(𝑡)𝑦𝑢(𝑦,𝑡)d𝑦𝑀𝑤,(4.1) with 𝑀𝑤 being the total amount of water. The integrals are evaluated numerically using the trapezoidal rule. If all the water is included, the superindex is dropped, so 𝑀𝑟 and 𝐺 are used; that is, these characteristics are linked with all the water in 𝑟(𝑟0(𝑡),𝑟0+𝐿+𝑅).

Experiment 1 (accuracy). The efficiency of the numerical method in solving the direct problem with sufficient accuracy is now demonstrated. The time evolution of the saturation and head profile is shown in 20 equidistant time sections for 𝑡(0,𝑇out) in Figures 2 and 3. For the same time sections 𝑠1, 𝑠2 and are presented in Figure 4. The corresponding values of global characteristics 𝑀𝑟, 𝐺, 𝑀𝑤(𝑝) in 10 time sections (𝑡𝑖=𝑖400,𝑖=1,,10, 𝑡11=𝑇out) are included in Table 1. Due to the numerical model satisfying the mass conservation principle, a global water mass 𝑀𝑤=6.1595 is found for the entire simulation time. The water amount of 0.1595 corresponds to an initial regularization of 𝑠1=0.2, 𝑠2=0.4 (zero values are singular for the ODE system); a smaller regularization can be used if needed. The entire computation only needs 10 seconds on a normal desktop PC. In Figure 3 only the nonpositive part of the head in the sample is shown, that is, only the region (𝑠1(𝑡),𝑠2(𝑡)). In the region (0,𝑠1(𝑡)), the head is a positive, convex parabola as given by formula (2.7).
Next, the centrifugation is continued over a time interval 𝑡(5016=𝑇out,9016) with outflow of the water. The corresponding saturation, saturation front, and expelled water level are drawn for 11 time sections in Figures 5 and 6. The global characteristics corresponding to these time sections are presented in Table 2. The following conclusions can be drawn. The global characteristics show significant change over the simulation time, indicating a good sensitivity. Figure 2 shows that the entire soil retention curve is needed to compute a solution during imbibition, and the use of a nonlocal boundary condition based on mass conservation allows simulation of the free outflow boundary, evident in Figure 5. A consequence of this is that for 𝑡>𝑇out, the flow is still governed partially by the unsaturated flow regime (𝑢0.5 at 𝑟=𝐿), making it useful to continue centrifugation in order to determine the soil parameters.
We now present how the soil parameters can be identified. For this, we focus on the parameters 𝛾 and 𝑛 in the van Genuchten model.

Table 1: Water in the sample, rotational momentum, and gravitational center, for Experiment 1.
Table 2: Rotational momentum and gravitational center.
Figure 4: Time evolution of the saturation front (full line), wetness front (dashed), and water level in the injection chamber (dash-dot, right scale) for 𝑡(0,5016=𝑇out).
Figure 5: Time evolution of the saturation 𝑡=5016+𝑖400,𝑖=0,,10.
Figure 6: Time evolution of the saturation front 𝑠1 (full line) and the expelled water 𝑤out (dashed line, right scale), 𝑡(5016,9016).

Experiment 2 (parameter identification). The data from Experiment 1 is used to construct an inverse problem: global measurements are extracted from this direct run and perturbed, and an inverse method (LM) is used to retrieve the original soil parameters. We focus on the parameters 𝛾, 𝑛, which in the direct simulation are. 𝛾=0.0189,𝑛=2.81.(4.2) During the time interval (0,𝑇out) the global characteristics that can be measured are [,𝑀𝑟,𝐺,𝑠2], while for 𝑡>𝑇out they are [,𝑀𝑟,𝐺,𝑤out]. Variations of these measurements will be considered to determine their influence. The wetness front 𝑠2 has been added here as it is a sharp interface that hence can be measured in several ways (not a requirement, see further). All measurements are collected in a vector we denote by 𝐃 the data. Practically, consider measurements at 10 uniformly distributed time sections with the time increment 𝑑𝑡=400 for 𝑡<𝑇out and 𝑑𝑡=200 for 𝑡>𝑇out. The end time is 𝑡=9016, while from the direct run we know 𝑇out=5016.
Before starting the inverse modeling, the data is perturbed by 5% or 10% by means of a random function. The resulting measurement is denoted by 𝐃5, 𝐃10, respectively, and is the approximation of the real measurement performed. The resulting soil parameters are indicated with the same subindex; for example, 𝛾5 is the resulting value for 𝛾 from 𝐃5. The LM method requires a starting parameter set, which we denote by 𝛾, 𝑛, and minimizes 𝐃𝑝𝐃(𝑖)2. Here, 𝐃𝑝 is the given perturbed measurement (e.g., 𝐃5) and 𝐃(𝑖) corresponds to 𝛾(𝑖), 𝑛(𝑖) in the 𝑖th LM iteration step.
For different measurements used we will depict the resulting saturation-capillary pressure and hydraulic permeability-capillary pressure curves as given by (2.2). The curves corresponding to 𝛾, 𝑛 will be given by a full line, 𝛾5, 𝑛5 with dash-dots, 𝛾10, 𝑛10 with dots, and 𝛾, 𝑛 with circles.
2A, Rotational Momentum and Gravitational Center.
Starting from 𝛾=0.0129, 𝑛=2, and using [,𝑀𝑟,𝐺,𝑤out], the following values are retrieved: 𝛾5=0.0183,𝑛5𝛾=2.7,10=0.0175,𝑛10=2.6,(4.3) corresponding to 𝐃5 and 𝐃10. The corresponding saturation-capillary pressure and hydraulic permeability-capillary pressure curves are given in Figure 7.
2B, Gravitational Center.
In the previous experiment the values of 𝑀𝑟 are dropped from the measurements, which gives the following values: 𝛾5=0.0188,𝑛5𝛾=2.79,10=0.0192,𝑛10=2.85.(4.4) Corresponding curves are given in Figure 8. We note that the reconstructed saturation-capillary pressure curve is better than when rotational moment is also considered. If in a similar experiment the values of 𝐺 are dropped instead of 𝑀𝑟, the restoration results are less good. This results corresponds to the observation that for the determined parameters the gravitational center is, in this specific case, more robust to the noise level than the rotational moment. Normally, adding more measurements improves the inverse method; however, this behavior was also seen in other experiments. Many more experiments would have to be done in order to make a general claim. This behavior is a weak point of the Levenberg-Marquardt method which includes the possibility of many local minima, see also Remark 4.1.
2C, Wetness Front.
We examine the influence of the evolution of the wetness front on the determination process. The measurements now consist of [,𝑠2,𝑤out]. We obtain 𝛾5=0.017,𝑛5𝛾=2.55,10=0.0166,𝑛10=2.49,(4.5) with corresponding curves in Figure 9. The result for the permeability curve is very good, but the saturation-capillary pressure curve is less good than in the other cases.
2D, Wetness Front and Gravitational Center.
We add now the 𝐺 measurements, giving data [,𝑠2,𝐺,𝑤out]. This gives 𝛾5=0.0175,𝑛5𝛾=2.6,10=0.017,𝑛10=2.54,(4.6) with corresponding curves in Figure 10. This result indicates a slight improvement as a consequence of adding the gravitational center measurements. However, this result is less good than when only the gravitational center is used, which again indicates it would be a good strategy to give the gravitational center a larger weight than the other measurements in the LM method.

We can conclude that the inverse determination works good, even with a perturbation of 10%. Using the levels in the water chambers combined with the gravitational center gives the best results. Errors on the rotational moment cause relatively large changes in the reconstructed saturation-capillary pressure curve. Changing the gravitational center measurements with those of the wetness front is an alternative but implies a penalty on the obtained values.

Figure 7: Data:[,𝑀𝑟,𝐺,𝑤out], (a) (𝑠/𝑝(𝑘)) curves, (b) (𝑝𝑒/𝑝(𝑘)) curves. Circles: start of inverse procedure, full line: correct result, dash-dots: result of inverse procedure with 5% noise, dots: result of inverse procedure with 10% noise.
Figure 8: Data:[,𝐺,𝑤out], (a) (𝑠/𝑝(𝑘)) curves, (b)(𝑝𝑒/𝑝(𝑘)) curves. Circles: start of inverse procedure, full line: correct result, dash-dots: result of inverse procedure with 5% noise, dots: result of inverse procedure with 10% noise.
Figure 9: Data:[,𝑠2,𝑤out], (a) (𝑠/𝑝(𝑘)) curves, (b) (𝑝𝑒/𝑝(𝑘)) curves. circles: start of inverse procedure, full line: correct result, dash-dots: result of inverse procedure with 5% noise, dots: result of inverse procedure with 10% noise.
Figure 10: Data:[,𝑠2,𝐺,𝑤out], (a) (𝑠/𝑝(𝑘)) curves, (b) (𝑝𝑒/𝑝(𝑘)) curves. Circles: start of inverse procedure, full line: correct result, dash-dots: result of inverse procedure with 5% noise, dots: result of inverse procedure with 10% noise.

Experiment 3 (saturated conductivity). We briefly show that also the saturated conductivity 𝐾𝑠 can easily be retrieved using the same model. For this, start with a saturated sample and water level (0)=3 in the injection chamber. Assume the measurements 𝐃 of the global characteristics [,𝑀𝑟,𝐺,𝑤out] in 10 uniformly distributed time section over (0,2000), and the perturbed data 𝐃5, 𝐃10, 𝐃15 as before. The starting value for the inverse model is 𝐾𝑠=105, while the correct one is 𝐾𝑠,=2.4105. We obtain with the LM method 𝐾𝑠,5=2.39105,𝐾𝑠,10=2.41105,𝐾𝑠,15=2.38105.(4.7) In the case where only 𝑀𝑟 measurements are used, we obtain 𝐾𝑠,5=2.35105,𝐾𝑠,10=2.44105,𝐾𝑠,15=2.33105.(4.8) Collecting only (𝑡) measurements, we obtain 𝐾𝑠,5=2.398105,𝐾𝑠,10=2.36105,𝐾𝑠,15=2.49105.(4.9) We can conclude that for the determination of 𝐾𝑠 all global characteristics lead to good identification, even in the presence of large perturbations.

Experiment 4 (drainage). We have determined the imbibition soil parameters in the second experiment. Now we present the setup to determine the drainage soil parameters. This is analogue with the well known outflow experiment, but we again augment the measurements with global characteristics. For drainage, a fully saturated sample is centrifugated over the time interval (0,10000𝑠). In the model this translates to a zero flux condition at the left boundary and controls outflow as before at the right boundary with outflow water collected in a water chamber. In a first experiment the measurements [𝑀𝑟,𝐺,𝑤out] from 20 uniformly distributed time sections are used. By LM iterations, starting from 𝛾=0.0129, 𝑛=2, we obtain 𝛾5=0.0176,𝑛5𝛾=2.68,10=0.0215,𝑛10=2.67.(4.10) Next, we investigate the influence of the rotational speed. Using now measurements from 20 uniformly distributed time sections in 𝑡(0,15000𝑠) from separate centrifugations with 3 different rotational speeds, we obtain the results from Table 3. The rotational speeds used are 𝜔=20, 30, and 40. Note that the average of these parameters (from 3 separate runs) gives a better estimate of the retention curve. We remark also that our experiments have shown that this average is better than the value obtained by LM iterations using the data collected in one centrifuge run where the rotational speed is changed after every 10000𝑠, taking consecutively Ω=20, 30, and 40.
As with the imbibition we now investigate if reducing the type of measurements still gives good results. Using the same measurements, but dropping the values of 𝑀𝑟, we obtain the results in Table 4. As in the imbibition case, the results are better. We remark that dropping the measurements of the gravitational center leads again to less good approximations of the parameters.
As a last numerical experiment, we apply consecutively two different rotational speeds 𝜔=20 and 𝜔=30, each over a period of 5000𝑠. Collecting as measurements 𝑀𝑟,𝐺,𝑤out over 10 uniformly distributed time sections along the interval 𝑡(0,10000𝑠), we obtain 𝛾5=0.0242,𝑛5𝛾=2.91,10=0.0192,𝑛10=2.9.(4.11) This result is worse than when only 𝜔=20 is used.
We can conclude from these experiments that the required soil parameters could be determined also in drainage mode with the proposed method.

Table 3: Estimated 𝛾, 𝑛 by LM method using all measurements for 3 rotational speeds.
Table 4: Estimated 𝛾, 𝑛 by LM method excluding rotational momentum measurements for 3 rotational speeds.

Experiment 5 (range of validity). In this experiment we examine the efficiency of the proposed algorithm on several different infiltration parameters. Consider the infiltration parameters 𝛾0=0.0189, 𝛾5=5𝛾0, 𝛾10=10𝛾0 and 𝐾𝑠,0=2.4.105, 𝐾𝑠,10=10𝐾𝑠,0, 𝐾𝑠,20=20𝐾𝑠,0, while 𝑛=2.81, Ω=30, and in the injection chamber we have initial water level 0=1.9. As before, there is a regularization so that at 𝑡=0 we have 𝑠1=0.2, 𝑠2=0.4. In Figures 11, 12, and 13, we draw the time evolution of the saturation distributions in 10 equidistant time points up to the time 𝑇𝑒, which is defined as the time at which the injection chamber is empty. This clearly shows that the algorithm can handle a wide range of parameter values.
Lastly, we perform some experiments to show that the regularization performed at the start for all computations up to now is not a limitation. For a practical experiment, one needs to start with a sample that is prepared and placed in the centrifuge. Hence, the sample will be in contact with water before the centrifugation starts. This allows to compute an initial profile of the saturation depending on this initial contact time. For simplicity we neglect gravitation here (it can be added to the algorithm in a similar fashion if needed). In Figures 14 and 15, the result is presented as two minutes of water infiltration into the sample from the water chamber in the absence of centrifugation (Ω=0). Here, a very small regularization at 𝑡=0 was performed corresponding with 𝑠1=0.01 and 𝑠2=0.02, which is small enough not to influence the final solution. From these figures it is clear that the contribution of infiltrated water (during 2 minutes) without centrifugation is significant, especially when the permeability is high. However, correct starting values for the full algorithm with Ω>0 can be computed without problems.

Figure 11: Time evolution of the saturation for 𝑡=𝑖𝑇𝑒/10, 𝑖=0,,10, (a) 𝛾0, 𝐾𝑠=𝐾𝑠,0, 𝑇𝑒=1516; (b) 𝛾0, 𝐾𝑠,10, 𝑇𝑒=167.
Figure 12: Time evolution of the saturation for 𝑡=𝑖𝑇𝑒/10, 𝑖=0,,10, (a) 𝛾5, 𝐾𝑠,10, 𝑇𝑒=187; (b) 𝛾10, 𝐾𝑠,10, 𝑇𝑒=190.
Figure 13: Time evolution of the saturation for 𝑡=𝑖𝑇𝑒/10, 𝑖=0,,10, with 𝛾5, 𝐾𝑠,20, 𝑇𝑒=98.
Figure 14: Time evolution of the saturation for Ω=0, 𝑡(0,120), (a) 𝛾0, 𝐾𝑠,0; (b) 𝛾5, 𝐾𝑠,10.
Figure 15: Time evolution of the saturation for Ω=0, 𝑡(0,120), 𝛾5, 𝐾𝑠,20.

Remark 4.1. As is typical, the inverse problem to determine 𝛾 and 𝑛 has (generally) many local minima. Applying the LM method, the obtained results (local minima) depend on the starting point chosen, and on several parameters inherent to the LM method. Therefore, one can run LM with different method parameters for the same starting point and then choose the “optimal’’ one defined by a minimal RMS. Next, one can create the set of such “optimal’’ solutions corresponding to different starting points and use a statistical analysis to create a reliable result. Sometimes it suffices to use an average of “local optimal’’ solutions. Between them we can choose those candidates for which RMS<0.4𝐸 (where 𝐸=𝐃𝐃102, a measure of the measurement error). We have applied this procedure in Experiment 4 with 4 local minima.

Remark 4.2. One can increase the reliability of the obtained 𝛾 and 𝑛 by using more measurements (extending the vector 𝐃 of measured characteristics). Also, in those cases where the imbibition/drainage hysteresis phenomenon can be neglected, the imbibition and drainage scenario above can be combined in one set of measurements, leading also to more measurements available. In those cases more complicated scenarios become also possible, like refilling the injection chamber after partial drainage of the sample.

5. Conclusions

A new and efficient numerical approximation for a water-air infiltration problem in porous media (in 1D) is developed, allowing for a free outflow boundary and resolution of the interfaces arising. It is shown with this model that it is sufficient to use only global characteristics, like rotational momentum, gravitational center, and inflow-outflow of water, to reconstruct the soil parameters in the governing mathematical model. The numerical experiments show that for practical realization, concentrating on measuring the gravitational center over the rotational momentum is advisable.

On the base of the developed numerical approximation, it is possible to use different centrifugation scenarios, which could be used for the realistic determination of soil parameters in imbibition or drainage mode.

It is also shown that to measure the wetness front in the sample is an alternative to measuring the rotational momentum and the gravitational center in the imbibition scenario. This front is very sharp during the centrifugation, allowing its precise determination with a detector (e.g., visually from the outside via a colorant and high-speed camera, or alternatively via X-rays to detect the water).

Based on the results of this numerical study, we propose to construct a centrifuge allowing to measure the gravitational center. Even if more detailed models are required (adding, e.g., consolidation, swelling, a 3D model with the Coriolis effect, …), the conclusion drawn here—several global characteristics are sufficient for determining the soil parameters—will remain valid.


The first and the third authors confirm financial support by the Slovak Research and Development Agency under Contract APVV-0184-10, APVV-0743-10, and VEGA-1/0502/09. The second author thanks the MaCKiE project (BOF-GOA 01GA0405) of the Ghent University which made the collaboration of the 3 authors possible.


  1. G. L. Hassler and E. Brunner, “Measurements of capillary pressure in small core samples,” AIME Transactions, vol. 160, pp. 114–123, 1945. View at Google Scholar
  2. J. R. Nimmo and K. A. Mello, “Centrifugal techniques for measuring saturated hydraulic conductivity,” Water Resources Research, vol. 27, no. 6, pp. 1263–1269, 1991. View at Publisher · View at Google Scholar · View at Scopus
  3. D. Constales and J. Kačur, “Determination of soil parameters via the solution of inverse problems in infiltration,” Computational Geosciences, vol. 5, no. 1, pp. 25–46, 2001. View at Publisher · View at Google Scholar · View at Scopus
  4. J. Šimůnek and J. R. Nimmo, “Estimating soil hydraulic parameters from transient flow experiments in a centrifuge using parameter optimization technique,” Water Resources Research, vol. 41, no. 4, pp. 1–9, 2005. View at Publisher · View at Google Scholar
  5. P. Forbes, “Discussion of techien chen, 1997, “An integral method to calculate water saturation in a centrifuge experiment to determine capillary pressure”,” Transport in Porous Media, vol. 31, no. 2, pp. 239–244, 1998. View at Google Scholar · View at Scopus
  6. S. Finsterle, T. O. Sonnenborg, and B. Faybishenko, “Inverse modeling of a multistep outflow experiment for determining hysteretic hydraulic properties,” in Proceedings TOUGH Workshop '98, pp. 250–256, Lawrence Berkeley National Laboratory, Berkeley, Calif, USA, 1998. View at Google Scholar
  7. I. S. Pop, C. J. van Duijn, J. Niessner, and S. M. Hassanizadeh, “Horizontal redistribution of fluids in a porous medium: the role of interfacial area in modeling hysteresis,” Advances in Water Resources, vol. 32, no. 3, pp. 383–390, 2009. View at Publisher · View at Google Scholar · View at Scopus
  8. T. Sakaki, D. M. O'Carroll, and T. H. Illangasekare, “Direct quantification of dynamic effects in capillary pressure for drainage-wetting cycles,” Vadose Zone Journal, vol. 9, no. 2, pp. 424–437, 2010. View at Publisher · View at Google Scholar · View at Scopus
  9. J. Kačur, B. Malengier, and P. Kišon, “A numerical model of transient unsaturated flow under centrifugation based on mass balance,” Transport in Porous Media, vol. 87, no. 3, pp. 793–813, 2011. View at Google Scholar
  10. M. T. van Genuchten, “A closed-form equation for predicting the hydraulic conductivity of unsaturated soils,” Soil Science Society of America Journal, vol. 44, no. 5, pp. 892–898, 1980. View at Google Scholar · View at Scopus
  11. C. E. Kees and C. T. Miller, “Higher order time integration methods for two-phase flow,” Advances in Water Resources, vol. 25, no. 2, pp. 159–177, 2002. View at Publisher · View at Google Scholar · View at Scopus