Numerical modelling of unsaturated-saturated flow under centrifugation with no outflow
A novel centrifuge set-up for the study of unsaturated flow characteristics in porous media is examined. In this set-up, simple boundary conditions can be used, but a free moving boundary between unsaturated-saturated flow arises. A precise and numerically efficient approximation is presented for the mathematical model based on Richards’ nonlinear and degenerate equation expressed in terms of effective saturation using the Van Genuchten-Mualem ansatz for the soil parameters in the unsaturated zone. Sensitivity of the measurable quantities (rotational moment, center of gravity and time period to achieve quasi steady state) on the soil parameters is investigated in several numerical experiments. They show that the set-up is suitable for the determination of the soil parameters via the solution of an inverse problem in an iterative way, excluding the saturated hydraulic conductivity. For this parameter, an existing simple centrifuge set-up is repeated and augmented with transient measurements.
To predict the flow and solute transport in soils one needs the soil hydraulic properties in terms of soil parameters. Once determined, these parameters can be used as input data in the governing mathematical model. For unsaturated flow, this model is given in terms of the saturation and the pressure head in Richards’ equation (see below), which is a nonlinear and degenerate parabolic equation. Furthermore, when part of the sample is saturated, free boundaries between the saturated zone and the partially saturated zone arise, as well as between the dry and the partially saturated zone. This is a major problem for many modeling approaches, leading to experimental set-ups that avoid the formation of these boundaries.
The soil retention and hydraulic permeability functions linking the saturation and pressure head for unsaturated flow are expressed using the Van Genuchten-Mualem ansatz by means of soil parameters. Measuring these soil parameters is usually time consuming and tedious, especially for low conductive porous media. Several set-ups based on centrifugation have been proposed to obtain a large acceleration of the processes involved, see [2, 4, 5, 8, 3] and citations therein. These techniques have several disadvantages. An approach which aims for a steady-state flow regime inside the centrifuge, [2, 4], requires expensive and/or complex apparatus, and obtains only a few water content versus conductivity measurements per run. Also transient set-ups based on keeping a top boundary at a fixed prescribed setting, , are expensive. The quasi-steady centrifuge (QSC) method, , is a much simpler technique (a slowly emptying reservoir at the top that is refilled when needed), but requires that the criterion for steadiness of flow through the sample is relaxed, leading to higher uncertainty in the obtained results.
The alternatives for determining conductivity with a steady-state flow, couple a transient flow with parameter estimation techniques, see eg. [3, 8]. In this way, the conductivity and retention curve can be determined inversely over a large saturation domain. These methods require experiments of some state variables which relate to the conductivity. One-step or multi-step outflow methods are common in column experiments. Pressures or suctions are applied at the top or the bottom of the sample, while the outflow, and optionally pressure head within the sample, are measured. These measurements are then used to estimate the hydraulic parameters. This technique is transferred to the centrifuge device in . However, measuring precise outflow rates is problematic in a centrifuge. Therefore, their method for determination of the soil parameters is based on transient water content measurements (coming from electrodes installed in the specimen). This can be used in an equilibrium analysis as well as in a transient analysis. In this setting, the specimen is fully saturated at the beginning, and saturation is maintained at the outflow via a short lip at the end of the container. Good results are obtained, but there remain some disadvantages to this technique: there are few measurements close to saturation, leading to a high error in the prediction of the conductivity close to saturation, the sample needs to be disturbed to introduce the electrodes, and there is a very long waiting time in order to achieve the equilibrium when the equilibrium analysis approach is used.
The main goal of this manuscript is to develop a precise numerical method enabling to determine the soil parameters (via solution of inverse problem) in a very simple way requiring very cheap measurements.
In this article we consider any partially saturated sample which is sealed at the right boundary (from the center of centrifuge) and has no inflow at the left boundary. The amount of infiltrated water is fixed and would be important initial information, but in the fact, we can drop this requirement and consider the initial water content as an additional parameter requiring determination. The only measurements required in our setting is the rotational momentum and the center of gravity of the sample in a set of equilibria corresponding to some predetermined rotational speeds. This can be replaced however by with transient measurements of these variables allowing to reduce the centrifuge operating time. These measurements are sufficient due to the fact that the saturation profile at the equilibria is not depending on the initial distribution of water in specimen, but only on its amount, which, when the right boundary of the sample is sealed, is identical in all equilibria.
To use this procedure, we have to face serious difficulties in the numerical modeling. The main one is that if the right side of the sample reaches effective saturation, an interface between partially saturated zone and saturated zone appears. This boundary is very difficult to control numerically, causing problems with the mass balance conservation which is very important in this set-up.
To reach the equilibrium is an infinite asymptotic process, but after some time (eg. 1-3 days for low conductive material) the change of the rotational momentum and of the center of gravity can no longer be measured. At that moment, the rotational speed is increased, and the system moves towards a new corresponding equilibrium. Note that even when equilibrium was not reached and a small error is present in the measurements of the rotational momentum and the center of gravity, this will not influence the error at the higher equilibrium level. This error depends only on the running time of centrifugation at the actual rotational speed. The differences between applied rotational speeds are chosen in such a way that that the differences in outputs (rotational momentum and center of gravity) are technically well distinguishable.
Next, the soil parameters and eventually the amount of originally infiltrated water, can be determined by minimizing a cost functional expressing the distance between the measured and the computed output, e.g., with the Levenberg-Marquard method. In this manuscript we only show the sensitivity of the outputs on the parameters, which is a good indication of how an inverse method will perform. In fact, the correct determination of soil parameters substantially depends on the sensitivity and the errors in practical measurements.
We can also use the transient measurements of the outputs, for example the period needed to cross from one equilibrium to another. This would give us additional information useful in the determination of the soil parameters.
The advantage of the above approach is that the full range of saturation values are present in the setup, while preventing outflow means equilibrium can be obtained faster. However, due to the set-up, it is clear that the water flows from the unsaturated zone to the saturated zone, with no flow occurring in the saturated zone. Indeed, we notice that the rotational momentum and center of gravity in the equilibria approach are not sufficiently sensitive on the “saturated hydraulic conductivity”. Therefore, we repeat in Section 4 also centrifugation with a fully saturated specimen which substantially increases sensitivity on this soil parameter, see . In this setup the sample remains saturated and the water table head in the top reservoir is measured over time and related to the saturated hydraulic conductivity. It is in other words the same as a column experiment where the dropping head at the top is measured instead of the outflow, and the head at the bottom is fixed instead of the head at the top.
In our numerical method we reduce the mathematical model to a system of ODE using method of lines (MOL), which has already been successfully applied to Richards’ equation in e.g. . Our main contribution is in correctly handling the moving free boundary. The obtained system can be solved with ODE solvers for stiff systems and stiff DAE systems (system of ODE and algebraic equations). The numerical experiments demonstrate that this solution method is numerically very efficient (precise and economical). Note that reaching equilibrium for a high rotational speed can require a week or even a month of centrifugation for low conductive materials, whereas the computational time is only a few seconds. Therefore this method is a good candidate for solving the inverse problem, as this requires many computations of the forward problem. The numerical method can be successfully applied in other centrifugation settings (concerning control of the inflow, or control of the outflow) as, e.g., in [8, 3].
The research performed in this article was done for a preliminary study for the feasibility of a new centrifuge device for low conductive materials that should be cheap to operate. The results presented here indicate that an approach where only rotational moment and gravity center is needed, is workable. This would remove the need to add experimental apparatus to the centrifuge measuring pressure heads, water contents or outflow rates during the operation. The authors intend to investigate a different set-up in a future article, in order to determine an optimal candidate, as well as collaborate to create a prototype centrifuge.
In Section 2 we present the mathematical model, giving specific attention to the movement of the free boundary. In Section 3 the numerical method based on the MOL approach is given, while in Section 4 the approach to determine the saturated hydraulic conductivity is explained. We finish in Section 5 with several numerical experiments showing the sensitivity of the output parameters on the soil parameters.
2 Mathematical model
We consider a one dimensional model for a partially saturated sample in the form of a tube. The tube starts (top or left boundary) at the distance from the center of the centrifuge and ends at the distance . The right boundary of the specimen is isolated. Flow in porous media under centrifugation is modeled by Darcy’s equation in the saturated region and by Richards’ equation in the unsaturated region (see, e.g., ,). So
in the saturated region, and
in the unsaturated region. Here, is the piesometric head, the saturation of the porous medium, the angular speed of rotation (in radians per second), the hydraulic conductivity in the saturated region, the gravitational constant and the function describes the hydraulic conductivity in the unsaturated region. Denote by the effective saturation, where is the volumetric water content at saturation and is the residual volumetric water content. We have , since . The soil hydraulic properties are represented by empirical expressions (see )
where , and are empirical soil parameters, is called the bubbling pressure .
It is possible to rewrite the flow in unsaturated form as
Equation (5) is strongly nonlinear and degenerate. We note that . Equilibria at the high rotational speed can be expected to have a fully saturated zone (supposing the initial amount of infiltrated water is sufficiently large), which appears at the right boundary and of which the front evolves to the left of the specimen (under non-decreasing rotational speed). We denote the position of this interface by . This saturated zone is governed by Darcy’s equation, but is unknown and time dependent. The time evolution of is difficult to compute. The dynamics of this region is linked with the (finite) interface flux
and based on a mass balance argument we can expect
Unfortunately, we cannot use this model for the determination of the time evolution of , since at it holds and . Consequently, .
If we transform Richards’ equation in terms of the head, we obtain
with , where is the hydraulic conductivity function,
and the specific moisture capacity function is given by
We can see that for . In Fig. 1 we present the graph of the functions and for , and parameter values , , .
As we can see also equation (7) degenerates at . This has to be taken into account when saturation becomes at the right boundary of specimen. After this moment, , the mathematical model must be changed to reflect the physical phenomenon. At the right hand side of the (isolated) specimen appears a saturated zone with an interface moving from the right boundary to the left. The flux at the interface is equal to , but also in this pressure-head form of Richards’ equation it is difficult to approximate correctly , which leads to a significant error in the mass balance.
Therefore to determine the interface , we will consider the algebraic equation
where is the amount of infiltrated water (which remains constant during the centrifugation). This condition reflects the global mass balance in the specimen and does not suffer from a flux approximation at .
Then mathematical model (7) only needs to be solved over the interval with right boundary condition for all . We approximate this mathematical model in the next section.
3 Numerical method
For the output parameters that will be measured (gravity center and rotational momentum), there is no need to model the head in the saturated zone, as we consider the compressibility of water to be negligible. The numerical approximation of mathematical model (7)-(8) results in a coupled system of a PDE and an algebraic equation. Moreover, the solution domain is a moving region, with unknown interface , which has to be determined.
We shift (7) to the domain and use the fixed domain transformation . This gives
Consider the space discretization
and , , and integrate (9) over for where , .
We denote by , , and approximate in the interval . We approximate
and similarly we approximate and denote it by . Let be the second order Lagrange polynomial crossing the points and . We use the abbreviation Then, our approximation of (9) (based on finite volume type approximation) at the point reads as follows
for . We add the corresponding equation at the point taking into account that the flux is zero there. In a similar way as in (10) (following the finite volume type of approximation) we obtain
At the point we have , so no additional equation is considered. We approximate the amount of water using the trapezoidal rule for the integration. Define
where . The last equation of this system is just (12). This system can be readily solved, e.g., by the solver “ode15s” in MATLAB or the “ida” solver of the Sundials package.
As is usual with these solvers, some regularization in (12) is needed as well as a tuning of the space discretization. Most important is to have a “good” starting point.
If the equilibria have the property , then no interface appears. It is then needed to set in the previous mathematical model and replace algebraic equation (12) by an ODE equation for which will be similar to (11). Successively increasing the rotational speed of the centrifuge increases the head at the right boundary. The model remains in the state where up to the point when , at which point the computation is automatically halted. The full model (10)-(12) is used onwards to compute the equilibrium states.
In the numerical equilibrium experiments (see Section 5, Exp. 5.5 and 5.7) it is observed, as expected, that the values of the rotational moment and the center of gravity are not very sensitive on the important parameter. Also the transient experiments where the time sections between different equilibria are measured, are not very sensitive. The saturated conductivity can only be determined from measurements of , that are accurate up to digits. Therefore, another method must be used for the determination of , which we present in the next section.
4 Centrifugation of saturated sample
For the determination of the saturated conductivity we propose to use the device put forward in , with the addition of allowing for transient measurements. We specifically use the ability to measure when a reservoir has completely drained out, combined with the measurements of the rotational moment.
We consider a saturated sample placed in the region of the centrifuge. A water reservoir is placed in front of the sample, and a head of is maintained at the bottom via a lip (alternatively one could place an empty reservoir into which the water flowing out from the sample is collected). The water from the front reservoir will flow through the sample and drain out to the right. We proceed the centrifugation up to the moment which is defined as the moment when the front reservoir is empty. By means of this setting the measurements of rotational momentum, resp. , will be very sensitive on the model parameter . The information of only the end time might even be sufficient to determine .
The mathematical model is very simple, since we have a saturated sample during the centrifugation - see (1), with time dependent Dirichlet boundary condition. The flux of water is given by
and is constant along all the sample. The head at the left boundary depends on the water level (in the reservoir) and is under centrifugation given by
The head at the outflow boundary is zero, i.e., . Hence, by simple integration of (1) in the saturated region we obtain the solution
Due to the mass balance argument we obtain the governing ODE for the time evolution of the water level:
with and . Solving this ODE we obtain the relation between and , whereas fully determines the change of the rotational moment over time. The numerical experiments for this set-up is presented in Sec. 5.8, below.
5 Numerical experiments
For all the experiments we use as data , , , , , , , , except where noted otherwise or where sequences are compared to investigate the sensitivity of the set-up on the parameters. If not mentioned otherwise, a uniformly distributed space discretization with grid points is used.
5.1 Independence of initial data
In the first numerical experiment we demonstrate that the equilibrium depends only on the rotational speed and the amount of infiltrated water. In Fig. 2 we show the time evolution of a set of unsaturated profiles over equidistant time sections in the interval . From this figure it is clear that at the end time the profile is very close to the equilibrium (of which the precise value is obtained in infinite time). Moreover, a long time before the end of the centrifugation at , the saturation profiles do not change visibly anymore.
5.2 Rotational speed
We investigate the sensitivity on the rotational speed . For this we monotonically increase in a stepwise way, from up to with increment , resulting in saturation profiles. At each rotational speed we reach the corresponding “good approximated” equilibrium. The centrifugation from one equilibrium to the next one takes .
In Fig. 4 we present the saturation profiles of the equilibria, together with one extra profile corresponding to the right boundary becoming fully saturated. The corresponding heads are presented in Fig. 5.
The extra profile is the sixth curve, which occurs during and has value . At that moment, the ODE solver is stopped, the system is changed to the system with a moving boundary, Eq. (7)-(8), and the solver is restarted with initial condition given by the last potential head (resp. saturation) profile.
The values at the equilibria of the rotational momentum , the center of gravity , and total water content are collected in Table 1. These are computed from the saturation profiles.
Two rotational moments are given, denotes the starting rotational momentum corresponding to the saturation distribution obtained from the previous rotational speed, and denotes the rotational momentum at the actual rotational speed in equilibrium. The formulas for , and at time are:
and are all evaluated numerically using the trapezoidal rule. Note that if then . The sensitivity of the measured quantities on the changing water content is very good. It remains to determine the contribution of the different parameters.
We note first however that careful analysis of the head profiles in Fig. 5 shows that these are not yet parabola at the end times of a run at fixed rotational speed, especially toward the left of the sample. This indicates that the equilibrium has not yet been reached. We investigate this further in the next Experiment.
5.3 Reaching equilibrium
To investigate the head profiles we start in this experiment from the equilibrium corresponding to and use a rotational speed . The centrifuge normally operates up to seconds. At that time, equilibrium for is almost reached. We compare values, the starting value, 9 increasing time steps (with ), the sensible end time step , and 2 extra time steps to investigate the very long time behavior. Fig. 6 shows the obtained saturation profiles, the corresponding evolution of the potential head is presented in Fig. 7.
As we can see the head profile at is still not a parabola in the left side of the profile (low head values). The reason for this is that the hydraulic permeability at low head is negligibly small, so it takes a very long time to reach the equilibrium. If the centrifugation is continued, also this part obtains the required parabolic shape associated with the equilibrium. Note however, that the other part of the head profile (for higher head values) is changing insignificantly. Therefore, we arrive at the conclusion that it makes sense to increase the rotational speed and not wait for these lower head values to stabilize. More so in light of Exp. 5.1.
5.4 Dependence on
In this experiment we demonstrate the sensitivity of and on the model parameter . We start with a constant saturation and apply now the rotational speed . The centrifuge is operated for . In Fig. 8 the obtained equilibrium profiles are depicted for successively ; ; ; ; ; ; and . In Fig. 9 the corresponding heads are given.
The resulting values for , and are given in Table 3, and indicate a good sensitivity.
5.5 Saturated conductivity
We now proceed with investigating the sensitivity of and on the saturated conductivity . We consider for ; ; ; ; ; and , start again from a uniform saturation, and apply this time a rotational speed of . In Table 4 we present the resulting values of , and at the obtained equilibrium.
The sensitivity of and on is negligible at the relative low pressure head that is generated for .
In Table 5 we present analogous results for the equilibrium obtained at rotational speed , starting from the equilibrium at , and considering a much larger spectrum for .
To reach the equilibrium from which Table 5 is extracted, a centrifugation time of is used for the values , whereas for values we have used an operating time of .
In Fig. 10 the equilibrium profiles of the saturation is given corresponding with the results from Table 5, while the corresponding heads are given in Fig. 11. From Fig. 11 one can see that for the lower values of the head a parabolic shape of the profile was still not reached if .
Overall, the results clearly indicate that another technique must be used to determine .
5.6 Dependence on
We now investigate the sensitivity of and on the last model parameter: . We again use a rotational speed of , starting from the equilibrium position at . As values for we consider with where increments of size are used. The values of and are listed in Table 6 and Table 7, respectively. The corresponding saturation and head profiles at time section are in Fig. 12 and Fig. 13, respectively.
5.7 Transient saturation determination
In experiment 5.5 we have shown the low sensitivity of on in an equilibrium analysis. One can assume this might improve in a transient analysis, as flow around the interface will be determined by . We investigate this now.
To assure the sample is in part fully saturated, high rotational speeds are used. At fixed we start from the equilibrium saturation at and let the centrifuge operate for (hours) with . After this time we change the rotational speed with the increment and again operate for the same duration. This is continued up to rotational speed , leading to end starting times , and five end times , .
The same procedure is applied for a series of values. Then we compare the values of and at the same time levels for different . In Table 8 we present the values of at time steps , (index corresponds to the rows). The columns correspond to the values with ; ; ; ; ; ; ; and .
So, in this table is computed s after each change of . The values of at time level , after each change of is presented in Table 9.
Although the changes in the tables are clear and consistent, they remain low. Tables 9 and 11 effectively correspond to the equilibrium measurements, whereas the comparison with Tables 8 and 10 gives an estimate of what is gained by a transient analysis. We can only conclude that for this set-up also a transient analysis to determine adds little improvements to the sensitivity. Note that for the parameters where an equilibrium analysis results in good sensitivity, a transient analysis is possible and can dramatically decrease the centrifuge operating time.
5.8 Saturated conductivity - 2
In this experiment we demonstrate that the determination of model parameter can be easily realized by means of very simple measurements obtained by centrifugation of a fully saturated sample with a water reservoir as presented in Section 4. In Fig. 14 and 15 the time evolution of the water level and of the rotational momentum is given. Note that in the computation of also a water collector placed at is considered. The parameter values considered are with ; ; ; ; and . The sensitivity of on will increase with higher values of .
The first and the third author confirm financial support by the Slovak Research and Development Agency under contract APVV-0351-07.
- M.C. Caputo and J.R Nimmo. Quasi-steady centrifuge method for unsaturated hydraulic properties. Water Resour. Res., 41:W11504, 2005.
- J.L. Conca and J.V. Wright. The ufa method for rapid, direct measurements of unsaturated transport properties in soil, sediment and rock. Aust. J. Soil Res., 36:291–315, 1998.
- D. Constales and J. Kačur. Determination of soil parameters via the solution of inverse problems in infiltration. Computational Geosciences, 5:25–46, 2004.
- R.J. Nimmo, K.C. Akstin, and K.A. Mello. Improved apparatus for measuring hydraulic conductivity at low water content. Soil Sci. Soc. Am. J., 56:1758–1761, 1992.
- R.J. Nimmo and K.A. Mello. Centrifugal techniques for measuring saturated hydraulic conductivity. Water Resour. Res., 27(6):1263–1269, 1991.
- M.D. Tocci and C.T. Kelley. Accurate and economical solution of the pressure-head form of richards’ equation by the method of lines. Advances in Water Resources, 20(1):1–14, 1997.
- M. Th. van Genuchten. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci.Soc.Am. J., 44:892–898, 1980.
- J. Šimůnek and J.R. Nimmo. Estimating soil hydraulic parameters from transient flow experiments in a centrifuge using parameter optimizetion technique. Water Resour. Res., 41, 2005.