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

Abstract

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, [6–8]. 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 𝑟=𝑟₀ from the center of the centrifuge and ends at the distance 𝑟=𝑟₀+𝐿, see Figure 1. A water chamber is placed at position 𝑟∈(𝑟₀−ℓ₀,𝑟₀)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:𝐾_𝑠𝜕_𝑟𝜕^𝑟ℎ−𝜔²𝜕𝑔𝑟=0,ℎ≥0,_𝑡𝜃=𝐾_𝑠𝜕_𝑟𝑘(𝜃)𝜕^𝑟ℎ−𝜔²𝑔𝑟,ℎ<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+(𝛾ℎ)^𝑛^𝑚𝑘(𝑢)=𝑢_1/2,ℎ∈(−∞,0),1−1−𝑢^1/𝑚^𝑚²,(2.2) where 𝑚=1−1/𝑛, 𝑛>1, and 𝛾=−(2^1/𝑚−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.3) where𝐷(𝑢)=−1(𝑛−1)𝛾𝜃_𝑠−𝜃_𝑟×𝑢^1/2−1/𝑚1−𝑢^1/𝑚^−𝑚1−1−𝑢^1/𝑚^𝑚².(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.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 𝑠₁(𝑡)∈(0,𝐿) the free boundary which separates the saturated and unsaturated regions in the sample, and denote by 𝑠₂(𝑡) the free boundary separating the partially saturated zone from the dry zone (wetness front), 𝑠₁(𝑡)<𝑠₂(𝑡)<𝐿. To obtain the solution for 𝑟∈(𝑟₀,𝑟₀+𝑠₁(𝑡)), the 2 integration parameters in (2.5) need to be determined from boundary conditions. We can easily determine thatℎ𝑟₀,𝑡=𝜔²2𝑔ℓ(𝑡)2𝑟⁰−ℓ(𝑡),ℎ𝑠₁(𝑡),𝑡=0,(2.6) since the piezometric head at 𝑟₀ is equal to the pressure of the water in the reservoir (contained in the interval (𝑟₀−ℓ(𝑡),𝑟₀)) and the piezometric head on the interface 𝑟₀+𝑠₁(𝑡) is zero. After elimination of the integration parameters we obtain the solution of (2.1) in the saturated region 𝑟∈(𝑟₀,𝑟₀+𝑠₁(𝑡)) in the formℎ(𝑟,𝑡)=𝜔²2𝑔𝑟−𝑟⁰²−𝑟−𝑟₀1𝑠₁ℓ2𝑟₀−ℓ+𝑠₁+ℓ2𝑟⁰−ℓ.(2.7) From this we can derive the flux 𝑞(𝑡) at the interface as𝑞(𝑡)=−𝐾_𝑠𝜕^𝑦ℎ𝑟₀+𝑦||_{𝑦=𝑠1(𝑡)}−𝜔²𝑔𝑟⁰=𝐾_𝑠+𝑠₁(𝑡)𝜔²2𝑔𝑠₁1(𝑡)2𝑟⁰𝑠₁(𝑡)+𝑠₁(𝑡)²+ℓ(𝑡)2𝑟₀−ℓ(𝑡),(2.8) where ℓ(𝑡) and 𝑠₁(𝑡) have to be determined. Since the flux of water along (𝑟₀,𝑟₀+𝑠₁(𝑡)) is constant (at fixed 𝑡) due to mass balance, it follows thatℓ̇(𝑡)=−𝑞(𝑡)∶=𝑓₁ℓ,𝑠₁,(2.9) where 𝑞 is obtained from (2.8). We proceed by constructing a model for 𝑠₁̇(𝑡). 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 𝑟∈(𝑟₀+𝑠₁(𝑡),𝑟₀+𝑠₂(𝑡)). We note that 𝑢(𝑠₂(𝑡),𝑡)=0 and 𝑢(𝑠₁(𝑡),𝑡)=1, forall𝑡. Then 𝑢 in the unsaturated region is governed by (2.1) with the moving domain of solution (𝑠₁(𝑡),𝑠₂(𝑡)) (we shift the original domain by 𝑟₀). We transform the moving domain to a fixed domain, using the transformation 𝑦=(𝑟−𝑠₁(𝑡))/(𝑠₂(𝑡)−𝑠₁(𝑡)). Denote 𝑢(𝑦,𝑡)∶=𝑢(𝑟,𝑡), and, for simplicity, in the following we drop the overline, writing shorthand 𝑢 for 𝑢(𝑦,𝑡). We have that𝜕_𝑡𝑢(𝑦,𝑡)=𝐾^𝑠𝑠₂₁(𝑡)²𝜕_𝑦𝐷(𝑢)𝜕^𝑦𝑢−𝑠₂₁(𝑡)𝑘(𝑢)𝜔²𝑔𝑟⁰+𝑠₁̇+𝑠₁(𝑡)+𝑦𝑠₂₁(𝑡)(𝑡)(1−𝑦)+𝑠₂̇(𝑡)𝑦1𝑠₂₁(𝑡)𝜕^𝑦𝑢,(2.10) where 𝑠₂₁(𝑡)∶=𝑠₂(𝑡)−𝑠₁(𝑡), and with boundary and initial conditions𝑢(0,𝑡)=1,𝑢(1,𝑡)=0;𝑢(𝑦,0)∶=𝑢⁰(𝑦).(2.11) Modeling the interface 𝑠₂(𝑡), 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.12) where 𝑝=1/2+1/𝑚. Therefore,𝑠₂̇(𝑡)=−𝐾^𝑠𝑚²𝑝𝜃_𝑠−𝜃_𝑟𝛾(𝑛−1)𝑠₂₁1(𝑡)𝜕^𝑦𝑢(1,𝑡)^𝑝,where𝑝=12+𝑚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 𝑞(𝑡)=𝑠₁̇(𝑡)+𝑞_int(𝑡), from which we obtain𝑠₁̇(𝑡)=𝑞(𝑡)−𝑞_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 (𝑠₂) 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 𝑠₂=𝐿, 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 𝑥∈(𝑠₁(𝑡),𝑠₂(𝑡)) which was transformed to a fixed domain 𝑦∈(0,1). Apply the space discretization0=𝑦₀<𝑦₁<⋯<𝑦_𝑖<⋯𝑦_𝑁=1,(3.1) and define 𝛼₀=0,𝛼_𝑖∶=𝑦_𝑖−𝑦_𝑖−1,𝑖=1,…,𝑁. This discretization corresponds to moving grid points 𝑥_𝑖(𝑡)=𝑟₀+𝑠₁(𝑡)+𝑦_𝑖𝑠₂₁(𝑡) 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²²¹×𝐷_𝑖+1/2𝜕⁺𝑢_𝑖−𝐷_𝑖−1/2𝜕⁻𝑢_𝑖−𝜔²𝑠₂₁×𝑘_𝑖+1/2𝑔+𝑠₁̇𝑟₀+𝑠₁+𝑦_𝑖+1/2𝑠₂₁−𝑘_𝑖−1/2𝑟₀+𝑠₁+𝑦_𝑖−1/2𝑠₂₁1−𝑦_𝑖+𝑠₂̇𝑦_𝑖𝑑ℒ𝑧;𝑦^𝑖𝑑𝑧||||_{𝑧=𝑦𝑖}ℓ̇(𝑡)=𝑓₁,(3.3)𝑠₁̇ℓ,𝑠₁,(3.4)𝑠₂̇(𝑡)=𝑓₁ℓ,𝑠₁−𝑓₂𝑠₁,𝑠₂,𝑢₁,𝑢₂,(3.5)(𝑡)=𝑓₃𝑠₁,𝑠₂,𝑢_𝑁−2,𝑢_𝑁−1,(3.6) where ℒ(𝑧;𝑦_𝑖) is the Lagrange polynomial of the second order crossing the points (𝑦_𝑖−1,𝑢_𝑖−1(𝑡)), (𝑦_𝑖,𝑢_𝑖(𝑡)), and(𝑦_𝑖+1,𝑢_𝑖+1(𝑡)), and function 𝑓₃ is given by𝑓₃𝑠₁,𝑠₂,𝑢_𝑁−2,𝑢_𝑁−1=−𝐾^𝑠𝑚²𝑝𝜃_𝑠−𝜃_𝑟𝛾(𝑛−1)𝑠₂₁1(𝑡)𝑑ℒ⁻𝑧;𝑦_𝑁𝑑𝑧||||_𝑧=1^𝑝,(3.7) where 𝑝=1/2+1/𝑚 and ℒ₋(𝑧;𝑦_𝑁) is the Lagrange polynomial crossing the points (𝑦_𝑁−2,𝑢_𝑁−2(𝑡)), (𝑦_𝑁−1,𝑢_𝑁−1(𝑡)), and (1,0). The function 𝑓₂ represents the numerical approximation of the flux 𝑞_int and reads as follows:𝑓₂𝑠₁,𝑠₂,𝑢₁,𝑢₂=𝐾^𝑠𝑠₂₁𝜃_𝑠−𝜃_𝑟𝑑ℒ⁺𝑧;𝑦₀𝑑𝑧||||_𝑧=0−𝜔²𝑔𝑠¹,(3.8) where ℒ₊(𝑧;𝑦₀) is the Lagrange polynomial crossing the points (0,0), (𝑦₁,ℎ₁), and (𝑦₂,ℎ₂), 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 𝑢₀=1, ℎ₀=0, and 𝐷(1)=+∞. We replace 𝑠₁̇, 𝑠₂̇ in (3.3) by 𝑓₁−𝑓₂ and 𝑓₃, respectively, to obtain the ODE system𝑤̇=𝐹(𝑡,𝑤),𝑤=𝑢₁,…,𝑢_𝑁−1,ℓ,𝑠₁,𝑠₂,(3.10) and solve it by an ODE solver for stiff systems. This system must be completed by the initial state 𝑤(0)=[𝑢₁(0),…,𝑢_𝑁−1(0),ℓ(0),𝑠₁(0),𝑠₂(0)]. Due to the first experimental scenario (infiltration into a dry region), all components of 𝑢 would be zero, and also 𝑠₁(0)=𝑠₂(0)=0 and ℓ(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 𝑠₂(0)>𝑠₁(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 𝑠₁ 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𝜃_𝑠−𝜃_𝑟𝑠₁(𝑡)+^𝑁_𝑖=0𝑐_𝑖𝑢_𝑖+𝑤^out(𝑡)=ℓ(𝑡)−ℓ(0),(3.11) where𝑐_𝑖=𝑠₂₁(𝑡)𝛼^𝑖+1+𝛼_𝑖𝑐2,𝑖=1,…,𝑁−1,₀=𝑠₂₁(𝑡)𝛼⁰2,𝑐^𝑁=𝑠₂₁(𝑡)𝛼⁰2.(3.12) This mathematical model can be used up to the time 𝑡=𝑇_out when 𝑠₂(𝑇_out)=𝐿. At that moment it is needed to replace 𝑠₂ by 𝑢_𝑁≈𝑢(1,𝑡) and put 𝑠₂≡𝐿; 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 by𝑤̇_out(𝑡)=−𝐾_𝑠𝐷𝑢^𝑁𝜕_𝑦ℎ𝑟₀+𝑦||_𝑦=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𝑧=ℎ₁,…,ℎ_𝑖0,𝑢_𝑖0+1,…,𝑢_𝑁−1,𝑠₁,𝑢_𝑁,ℓ,𝑤_out.(3.15) In this, the index 𝑖₀ 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 𝑖₀=𝑁/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 𝑢₀. The governing ODE system is obtained from (3.3) with 𝑖=0,…,𝑁−1, where 𝑠₁=0, 𝑠₂=𝐿, and 𝑠₂₁=𝐿, ℓ=0. An ODE is added for 𝑢̇₀ similarly as for 𝑢̇_𝑁 before. The outflow 𝑤_out is governed by (3.13), where 𝑟₀+𝑠₁(𝑡) is replaced with 𝑟₀+𝐿. The system is completed with the mass balance equation (3.11), where 𝑠₁=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+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 byℓ̇(𝑡)=𝐾_𝑠𝜔²2𝑔𝐿𝐿²−ℓ(𝑡)²+2𝑟₀(𝐿+ℓ(𝑡)),(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: 𝑟₀=30, 𝐿=10, 𝜔=20, 𝐾_𝑠=2.410⁻⁵, 𝜃_𝑟=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 𝑑₁=𝐿/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 ℓ₀=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𝑀^𝑟(𝑝)(𝑡)=𝑠²¹(𝑡)𝜔²2¹⁰𝑟₀𝑀+𝑠₁+𝑠₂₁(𝑡)𝑦²𝑢(𝑦,𝑡)d𝑦,_𝑤(𝑝)𝐺(𝑡)=𝑠₂₁(𝑡)¹⁰𝑢(𝑦,𝑡)d𝑦,_(𝑝)(𝑡)=𝑠₂₁(𝑡)¹⁰𝑟₀+𝑠₁+𝑠₂₁(𝑡)𝑦𝑢(𝑦,𝑡)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 𝑟∈(𝑟₀−ℓ(𝑡),𝑟₀+𝐿+𝑅).

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 𝑠₁, 𝑠₂ and ℓ are presented in Figure 4. The corresponding values of global characteristics 𝑀_𝑟, 𝐺, 𝑀^𝑤(𝑝) in 10 time sections (𝑡_𝑖=𝑖⋅400,𝑖=1,…,10, 𝑡₁₁=𝑇_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 𝑠₁=0.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 (𝑠₁(𝑡),𝑠₂(𝑡)). In the region (0,𝑠₁(𝑡)), 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 [ℓ,𝑀_𝑟,𝐺,𝑠₂], while for 𝑡>𝑇_out they are [ℓ,𝑀_𝑟,𝐺,𝑤_out]. Variations of these measurements will be considered to determine their influence. The wetness front 𝑠₂ 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 𝐃₅, 𝐃₁₀, respectively, and is the approximation of the real measurement performed. The resulting soil parameters are indicated with the same subindex; for example, 𝛾₅ is the resulting value for 𝛾 from 𝐃₅. The LM method requires a starting parameter set, which we denote by 𝛾, 𝑛, and minimizes ‖𝐃_𝑝−𝐃^(𝑖)‖². Here, 𝐃_𝑝 is the given perturbed measurement (e.g., 𝐃₅) 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, 𝛾₅, 𝑛₅ with dash-dots, 𝛾₁₀, 𝑛₁₀ with dots, and 𝛾, 𝑛 with circles.
2A, Rotational Momentum and Gravitational Center.
Starting from 𝛾=−0.0129, 𝑛=2, and using [ℓ,𝑀_𝑟,𝐺,𝑤_out], the following values are retrieved: 𝛾𝛾₅=−0.0183,𝑛₅=2.7,₁₀=−0.0175,𝑛₁₀=2.6,(4.3) corresponding to 𝐃₅ and 𝐃₁₀. 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: 𝛾𝛾₅=−0.0188,𝑛₅=2.79,₁₀=−0.0192,𝑛₁₀=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 [ℓ,𝑠₂,𝑤_out]. We obtain 𝛾𝛾₅=−0.017,𝑛₅=2.55,₁₀=−0.0166,𝑛₁₀=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 [ℓ,𝑠₂,𝐺,𝑤_out]. This gives 𝛾𝛾₅=−0.0175,𝑛₅=2.6,₁₀=−0.017,𝑛₁₀=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.

(a)

(b)

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.

(a)

(b)

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.

(a)

(b)

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.

(a)

(b)

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 𝐃₅, 𝐃₁₀, 𝐃₁₅ as before. The starting value for the inverse model is 𝐾_𝑠=10⁻⁵, while the correct one is 𝐾_𝑠,⋆=2.410⁻⁵. We obtain with the LM method 𝐾_𝑠,5=2.39𝐾10⁻⁵_𝑠,15,𝐾_𝑠,10=2.4110⁻⁵,=2.3810⁻⁵.(4.7) In the case where only 𝑀_𝑟 measurements are used, we obtain 𝐾_𝑠,5=2.35𝐾10⁻⁵_𝑠,15,𝐾_𝑠,10=2.4410⁻⁵,=2.3310⁻⁵.(4.8) Collecting only ℓ(𝑡) measurements, we obtain 𝐾_𝑠,5=2.398𝐾10⁻⁵_𝑠,15,𝐾_𝑠,10=2.3610⁻⁵,=2.4910⁻⁵.(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 𝛾𝛾₅=−0.0176,𝑛₅=2.68,₁₀=−0.0215,𝑛₁₀=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 𝛾𝛾₅=−0.0242,𝑛₅=2.91,₁₀=−0.0192,𝑛₁₀=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.0189, 𝛾₅=5𝛾₀, 𝛾₁₀=10𝛾₀ and 𝐾_𝑠,0=2.4.10⁻⁵, 𝐾_𝑠,10=10𝐾_𝑠,0, 𝐾_𝑠,20=20𝐾_𝑠,0, while 𝑛=2.81, Ω=30, and in the injection chamber we have initial water level ℓ₀=1.9. As before, there is a regularization so that at 𝑡=0 we have 𝑠₁=0.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 𝑠₁=0.01 and 𝑠₂=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.

(a)

(b)

Figure 11 

Time evolution of the saturation for 𝑡=𝑖𝑇𝑒/10, 𝑖=0,…,10, (a) 𝛾0, 𝐾𝑠=𝐾𝑠,0, 𝑇𝑒=1516; (b) 𝛾0, 𝐾𝑠,10, 𝑇𝑒=167.

(a)

(b)

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.

(a)

(b)

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 𝐸=‖𝐃_⋆−𝐃₁₀‖², 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.

Acknowledgments

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.