Numerical modelling of unsaturated-saturated flow under centrifugation with no outflow

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, [3], are expensive. The quasi-steady centrifuge (QSC) method, [1], 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 [8]. 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 [5]. 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. [6]. 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.

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 $r=r_{0}$ from the center of the centrifuge and ends at the distance $r=r_{0}+L$. 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., [8],[3]). So

in the saturated region, and

in the unsaturated region. Here, $h$ is the piesometric head, $\theta$ the saturation of the porous medium, $\omega$ the angular speed of rotation (in radians per second), $K_{s}$ the hydraulic conductivity in the saturated region, $g$ the gravitational constant and the function $k(\theta)$ describes the hydraulic conductivity in the unsaturated region. Denote by $u=\frac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}}$ the effective saturation, where $\theta_{s}$ is the volumetric water content at saturation and $\theta_{r}$ is the residual volumetric water content. We have $u\in(0,1)$, since $\theta\in(\theta_{s},\theta_{r})$. The soil hydraulic properties are represented by empirical expressions (see [7])

where $m=1-1/n$, $n>1$ and $\gamma=-(2^{1/m}-1)^{1-m}/h_{b}$ are empirical soil parameters, $h_{b}$ is called the bubbling pressure .

It is possible to rewrite the flow in unsaturated form as

where

Equation (5) is strongly nonlinear and degenerate. We note that $D(0)=0,D(1)=\infty$. 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 $s(t)$. This saturated zone is governed by Darcy’s equation, but $s(t)$ is unknown and time dependent. The time evolution of $s(t)$ is difficult to compute. The dynamics of this region is linked with the (finite) interface flux $q_{i}$

and based on a mass balance argument we can expect

Unfortunately, we cannot use this model for the determination of the time evolution of $s(t)$, since at $r=s(t)$ it holds $u=1$ and $D(1)=\infty$. Consequently, $\partial_{r}u|_{r=s(t)}=0$.

If we transform Richards’ equation in terms of the head, we obtain

with $k_{0}=\frac{K_{s}}{\theta_{s}-\theta_{r}}$, where $k_{0}\bar{k}(h)$ is the hydraulic conductivity function,

and the specific moisture capacity function $d_{s}(h)=\mathrm{d}u/\mathrm{d}h$ is given by

We can see that $\bar{k}(h)\to 1$ for $h\to 0$. In Fig. 1 we present the graph of the functions $\bar{k}(h)$ and $100\,d_{s}(h)$ for $h\in(-200,0)$, and parameter values $K_{s}=2.4\,10^{-5}$, $n=2.81$, $\gamma=-0.0189$.

As we can see also equation (7) degenerates at $h=0$. This has to be taken into account when saturation becomes $1$ at the right boundary of specimen. After this moment, $t=t_{1}$, 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 $s(t)$ moving from the right boundary to the left. The flux at the interface $s(t)$ is equal to $-\dot{s}(t)$, but also in this pressure-head form of Richards’ equation it is difficult to approximate correctly $\partial_{t}h|_{x=s(t)}$, which leads to a significant error in the mass balance.

Therefore to determine the interface $s(t)$, we will consider the algebraic equation

where $M_{w}$ 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 $r=r_{0}+s(t)$.

Then mathematical model (7) only needs to be solved over the interval $r\in(r_{0},r_{0}+s(t))$ with right boundary condition $h(r_{0}+s(t))=0$ for all $t$. We approximate this mathematical model in the next section.

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 $s(t)$, which has to be determined.

We shift (7) to the domain $r\in(0,s(t))$ and use the fixed domain transformation $y=\frac{r}{s(t)}$. This gives

Consider the space discretization

and $\alpha_{0}=0$, $\alpha_{i}:=y_{i}-y_{i-1}$, $i=1,\ldots,N$ and integrate (9) over $I_{i}:=(y_{i-1/2},y_{i+1/2})$ for $i=1,\ldots,N-1$ where $y_{i-1/2}:=(y_{i}+y_{i-1})/2$, $y_{i+1/2}:=(y_{i}+y_{i+1})/2$.

We denote by $h_{i}(t)\approx h(y_{i},t)$, $\forall i=1,\ldots,N-1$, and approximate $\mathrm{d}_{t}h(y,t)\approx\dot{h}_{i}(t)$ in the interval $I_{i}$. We approximate

and similarly we approximate $\partial_{y}h|_{y=y_{i-1/2}}$ and denote it by $:\partial^{-}h_{i}$. Let ${\cal L}(z;y_{i})$ be the second order Lagrange polynomial crossing the points $(y_{i-1},h_{i-1}),(y_{i},h_{i})$ and $(y_{i+1},h_{i+1})$. We use the abbreviation $k_{i+1/2}:=\bar{k}(h_{y_{i+1/2}})$ Then, our approximation of (9) (based on finite volume type approximation) at the point $y=y_{i}$ reads as follows

for $i=1,...,N-1$. We add the corresponding equation at the point $y_{0}$ 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 $y_{N}=1$ we have $h_{N}(t)=0$, so no additional equation is considered. We approximate the amount of water $M_{w}$ using the trapezoidal rule for the integration. Define

where

Then, system (10)-(11) will be completed by the algebraic equation

This algebraic equation is used instead of an ODE equation that models $\dot{s}(t)$. System (10)-(12) is degenerate and is of the form

where $z=(h_{0},h_{1},...,h_{N-1},s)$. 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 $h_{N}<0$, then no interface appears. It is then needed to set $s(t):=L$ in the previous mathematical model and replace algebraic equation (12) by an ODE equation for $\dot{h}_{N}$ 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 $s(t):=L$ up to the point when $h(N)=0$, 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 $M_{r}$ and the center of gravity $G_{c}$ are not very sensitive on the important $K_{s}$ parameter. Also the transient experiments where the time sections between different equilibria are measured, are not very sensitive. The saturated conductivity $K_{s}$ can only be determined from measurements of $M_{r}$, $G_{c}$ that are accurate up to $3$ digits. Therefore, another method must be used for the determination of $K_{s}$, which we present in the next section.

For the determination of the saturated conductivity we propose to use the device put forward in [5], 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 $(r_{0},r_{0}+L)$ of the centrifuge. A water reservoir is placed in front of the sample, $r\in(r_{0}-l_{0},r_{0})$ and a head of $0$ is maintained at the bottom $r_{0}+L$ via a lip (alternatively one could place an empty reservoir $r\in(r_{0}+L,r_{0}+L+d)$ 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 $T_{e}$ which is defined as the moment when the front reservoir is empty. By means of this setting the measurements of rotational momentum, resp. $T_{e}$, will be very sensitive on the model parameter $K_{s}$. The information of only the end time $T_{e}$ might even be sufficient to determine $K_{s}$.

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 $r=r_{0}$ depends on the water level $\ell(t)$ (in the reservoir) and is under centrifugation given by

The head at the outflow boundary is zero, i.e., $h(r_{0}+L)=0$. Hence, by simple integration of (1) in the saturated region we obtain the solution

where $x\in(0,L)$. Inserting (15) in (14) we obtain

Due to the mass balance argument we obtain the governing ODE for the time evolution of the water level:

with $\ell(0)=l_{0}$ and $\ell(T_{e})=0$. Solving this ODE we obtain the relation between $T_{e}$ and $K_{s}$, whereas $\ell(t)$ fully determines the change of the rotational moment $M_{r}(t)$ over time. The numerical experiments for this set-up is presented in Sec. 5.8, below.

For all the experiments we use as data $r_{0}=10$, $L=10$, $\omega=30$, $K_{s}=2.4\,10^{-5}$, $\theta_{r}=0.02$, $\theta_{s}=0.4$, $\gamma=-0.0189$, $n=2.81$, 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 $N=40$ grid points is used.

The first and the third author confirm financial support by the Slovak Research and Development Agency under contract APVV-0351-07.