Mathematical Problems in EngineeringVolume 2011 (2011), Article ID 163020, 23 pagesdoi:10.1155/2011/163020

Research Article

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

Jozef Kačur,¹ Benny Malengier,² and Pavol Kišon¹

¹Faculty of Mathematics, Physics and Informatics, Comenius University Bratislava, Mlynská dolina, 84248 Bratislava, Slovakia²Department 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.

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 , and outflow chamber to the right. Full saturation up to and dry zone to the right of .

Figure 2: Time evolution of the saturation, . Left curve at , right curve for .

Saturated flow (piezometric head ) in porous media under centrifugation is modeled by Darcy’s equation, while the unsaturated part () is governed by Richards’ equation as follows: 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 , since . We consider soil hydraulic properties which were proposed in [10] where , , and are empirical soil parameters and is the bubbling pressure. Hence, flow in an unsaturated region can be rewritten in terms of saturation as where Equation (2.3) is strongly nonlinear and degenerate (, ). 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 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 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 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 From this we can derive the flux at the interface as where and have to be determined. Since the flux of water along is constant (at fixed ) due to mass balance, it follows that 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 entering the unsaturated subregion. The unsaturated region is placed between two free boundaries . We note that and , . 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 where , and with boundary and initial conditions Modeling the interface , we follow [3] where problems similar to (2.10) have been studied. We compute first the order of degeneracy of : where . Therefore,

Now, the solution of (2.10)–(2.13) defines a flux which is the flux of the solution at . Finally, we close our system with the mass balance condition , from which we obtain 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 (), as in this scenario we want to avoid drainage (which starts from the left when ). If however at some we have 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 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 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 . Apply the space discretization and define . 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 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 . Denote and in the interval , . Approximate where , and similarly approximate and denote it by . Also denote and (similarly for ). Then, the approximation of (2.10) at the point , , reads as follows: where is the Lagrange polynomial of the second order crossing the points , , and, and function is given by where and is the Lagrange polynomial crossing the points , , and . The function represents the numerical approximation of the flux and reads as follows: where is the Lagrange polynomial crossing the points , , and , with We note that the flux at must be expressed in terms of head (instead of saturation) since at we have , , and . We replace , in (3.3) by and , respectively, to obtain the ODE system and solve it by an ODE solver for stiff systems. This system must be completed by the initial state . Due to the first experimental scenario (infiltration into a dry region), all components of would be zero, and also and . 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 , with in between.

Figure 3: Time evolution of the head, . Left curve at , right curve for .

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 where This mathematical model can be used up to the time when . At that moment it is needed to replace by and put ; that is, we drop (3.6) and add the ODE for to (3.3) (created on the same argumentation as for ). Moreover, the ODE system is completed by linking the outflow flux and . Now, the approximation consists of ODE and the algebraic equation (3.11). This system can be written in the form where 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 . 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 , where , , and , . An ODE is added for similarly as for before. The outflow is governed by (3.13), where is replaced with . The system is completed with the mass balance equation (3.11), where and 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

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 and we can proceed centrifugation up to .

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: , , , , , , , and . For the space discretization grid points are used with geometrical distribution: a first space interval with width , and then () with . 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 . At the start it is assumed that no water is present in the outflow chamber (which remains the case up to ). In the case of a filled injection chamber, initial water level 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 (as the surface is 1). The partially saturated zone was transformed to the domain . The contribution of the mobile water in the partially saturated zone to the global characteristics is denoted by superindex and can be calculated from 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 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 (, ) are included in Table 1. Due to the numerical model satisfying the mass conservation principle, a global water mass is found for the entire simulation time. The water amount of 0.1595 corresponds to an initial regularization of , (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 , the head is a positive, convex parabola as given by formula (2.7).Next, the centrifugation is continued over a time interval 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 , the flow is still governed partially by the unsaturated flow regime ( 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.

1
0.9100
0.9665
27.9404

2
1.3066
1.0031
28.4539

3
1.6256
1.0354
28.8971

4
1.9044
1.0658
29.3067

5
2.1576
1.0952
29.6965

6
2.3928
1.1241
30.0739

7
2.6144
1.1527
30.4430

8
2.8253
1.1813
30.8065

9
3.0277
1.2099
31.1662

10
3.2229
1.2386
31.5236

11
3.6931
1.3125
32.4262

Table 2: Rotational momentum and gravitational center.

1
1.3446
33.4371

2
1.3807
33.8598

3
1.4155
34.2633

4
1.4491
34.6520

5
1.4813
35.0254

6
1.5123
35.3844

7
1.5420
35.7296

8
1.5706
36.0617

9
1.5981
36.3808

10
1.6244
36.6867

11
1.6855
37.3956

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 .

Figure 5: Time evolution of the saturation .

Figure 6: Time evolution of the saturation front (full line) and the expelled water (dashed line, right scale), .

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. During the time interval the global characteristics that can be measured are , while for they are . 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 for and for . The end time is , while from the direct run we know .Before starting the inverse modeling, the data is perturbed by or 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 , , and using , the following values are retrieved: 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: 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 . We obtain 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 . This gives 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 . 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:, (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:, (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:, (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:, (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 in the injection chamber. Assume the measurements of the global characteristics in 10 uniformly distributed time section over , and the perturbed data , , as before. The starting value for the inverse model is , while the correct one is . We obtain with the LM method In the case where only measurements are used, we obtain Collecting only measurements, we obtain 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 . 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 from 20 uniformly distributed time sections are used. By LM iterations, starting from , , we obtain Next, we investigate the influence of the rotational speed. Using now measurements from 20 uniformly distributed time sections in from separate centrifugations with 3 different rotational speeds, we obtain the results from Table 3. The rotational speeds used are , 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 , taking consecutively , 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 and , each over a period of . Collecting as measurements over 10 uniformly distributed time sections along the interval , we obtain This result is worse than when only 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.

20

2.82

2.6

30

2.79

3

40

2.99

3

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

20

2.64

2.37

30

2.98

2.35

40

2.78

2.7

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 , , and , , , while , , and in the injection chamber we have initial water level . As before, there is a regularization so that at we have , . 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 (). Here, a very small regularization at was performed corresponding with and , 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 can be computed without problems.

(a)

(b)

Figure 11: Time evolution of the saturation for , , (a) , , ; (b) , , .

(a)

(b)

Figure 12: Time evolution of the saturation for , , (a) , , ; (b) , , .

Figure 13: Time evolution of the saturation for , , with , , .

(a)

(b)

Figure 14: Time evolution of the saturation for , , (a) , ; (b) , .

Figure 15: Time evolution of the saturation for , , , .

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 (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.

References

G. L. Hassler and E. Brunner, “Measurements of capillary pressure in small core samples,” AIME Transactions, vol. 160, pp. 114–123, 1945.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.