Inverse determination of saturated and relative permeability with a bench-scale centrifuge

Abstract

It is well established that the centrifuge can be an important tool to determine permeabilities of porous materials. Nevertheless, up to this date, for geotechnical applications, they are not used outside research laboratories. We attribute this to the cost of a centrifuge, the apparent requirement of measurements in the sample and the more complicated handling of samples for use in the centrifuge. As a first step towards an affordable centrifuge, we present an adaptation to a common bench-scale centrifuge, which allows to obtain measurements in a fraction of the time needed by standard infiltration methods, especially for clayey soils. The permeabilities are determined by measuring global transient data only: water level in an inflow or outflow chamber and gravitational centre of the sample or the rotational momentum of the sample. These are all measured outside of the centrifuge in a start-stop regime. A method of lines approach to solve the Richards’ equation is used, combined with interface tracking, which was previously shown to be an accurate method. Permeabilities are assumed to behave according to the van Genuchten–Mualem formula. We present the experimental set-up, accuracy of the measurements and the numerical solution method. For saturated flow, we compare with standard tests, while for unsaturated flow we compare with a pressure plate result (ASTM D2325). The results indicate the strengths and limitations of this approach, but overall we can conclude that measuring the global characteristics during rotation should provide accurate permeabilities.

Keywords:

• centrifuge
• inverse problem
• soil parameters
• unsaturated flow
• van Genuchten–Mualem ansatz

Classification:

• 65M20
• 65M32

Introduction

When predicting the performance of infiltration zones or clay covers, it is important to also determine the behaviour under unsaturated conditions. For example, compacted clay covers and linings of waste-management facilities are often in unsaturated conditions, and the same is true for infiltration zones after a long dry period. The clay covers and linings have the purpose to contain polluted liquids that might spread into the surrounding soil and groundwater. To predict the flow, and also the solute transport, through soils requires the determination of soil hydraulic properties. This means the saturated conductivity as well as the unsaturated conductivity mainly depends on the water content, dry density and degree of saturation of the soil. Many techniques exist to determine the unsaturated conductivity, see e.g. [1–4].

Various researchers have tried to evaluate unsaturated hydraulic conductivity (relative permeability) of soils by conducting either laboratory experiments or in situ studies. These methods can be classified as direct methods or indirect methods. The direct methods are quite tedious and time consuming and require expensive experimental set-ups.[5] Indirect methods, which employ volumetric properties of the soil and the water retention curve (also called soil-water characteristic curve, SWCC), are frequently adopted. Integration along the water retention curve provides a measure of the quantity of water in the soil which can then be used to estimate the soil hydraulic conductivity.[4, 6, 7] Although some field and laboratory methods are available to measure the relationship between water permeability and suction of unsaturated soils,[8] it is extremely difficult and time-consuming to measure the relationship accurately, especially for low conductive materials.

The determination of the permeability can be accelerated by using centrifugation. The measurement of capillary pressure curves with a centrifuge was initiated in [9]. 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 [10–12] (see also citations there).

Originally, for geotechnical uses, the majority of applications of centrifugation were 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. One approach applies the series of equilibria and the amount of expelled water to obtain the retention curve, see overview in [13]. This approach is relatively slow however. Recently, in [14], such a procedure was used to compare centrifuge results with a standard pressure plate method. A start-stop procedure was used to measure the water content by weight outside the centrifuge when equilibrium was obtained (120 min of centrifugation). Excellent agreement was found, while the testing time was reduced from 2 to 3 months with the pressure plate method to 4 days with the centrifuge (including oven drying).

A faster approach requires transient data, for example, from inside the sample at one or more points. This approach is similar to the common multi-step outflow experiments to determine the soil retention curve. For example,[15] shows that such an approach combined with an inverse method allows to determine hysteretic hydraulic properties of the soil retention curve. In [16], the distribution of the saturation in equilibria (linked with the corresponding rotational speeds) is measured via electrical signals from electrodes installed in the sample, and transient data used to find the hydraulic parameters. In general, the agreement is that for one step centrifuge experiments, one can use outflow water measurements combined with some internal measurements to obtain the soil retention curve.[17] Otherwise, the parameters are not unique. It is also well established that multi-rate experiments can decrease the uncertainty in the estimated hydraulic functions.[12] Using outflow water measurements seems to overestimate the soil retention curve, though this was not found in [14]. This means that the soil is more permeable than what is determined via the centrifuge method.

Although the benefits of using the centrifuge are well known, they are not yet used by geotechnical laboratories to determine permeabilities. We attribute this to three factors. First, the used centrifuges are often expensive. This is because they are large (arms $>$1 m) or because large modification is needed. E.g. [17] use image recognition in the outflow cup. For actual use, the price of a immediately usable centrifuge should not be above 10.000 Euro. A second reason is that operation is more complicated than using the classical pressure plate method. Easy-to-understand instructions are needed, and the user needs to have the feeling the centrifuge will lead to clear results. When transient measurements are used, operators need to rely on a black box numerical method for the results, which is not an approach they are accustomed to. Hence, strong proof of the centrifuge method works accurately is needed. As a third reason, the need for internal measurements is a draw-back. When internally tensiometers must be used, sample preparation is not straightforward.

In a series of articles,[18, 19] it was shown mathematically that transient measurements of only global characteristics in a centrifuge would be sufficient to determine the soil parameters governing the conductivity. This could greatly decrease the complexity of the centrifuge used, reducing the handling complexity and the centrifuge cost. In this paper, we show experimentally that the mathematical claims of these articles are true, but that further adaptations to the centrifuge are needed. Specifically, a cheap centrifuge is used, with minimal adaptations to perform the experimental work. A start-stop approach is used to investigate the claims, measuring in a transient way outflow water and gravitational centre of the sample. In order for the start-stop scheme to work, the model was adapted to allow for a hydrostatic pressure boundary like in [17].

Innovative is that a table top centrifuge is used, and a transient mathematical model is applied to determine inversely the soil parameters arising in the governing equation (Richards’ equation) and the van Genuchten–Mualem formula. In order to check the precision of the centrifuge, the saturated permeability is first determined for filters and samples. Next, the relative permeability of the soil samples is determined. We conclude this paper with recommendations on how to further adapt the centrifuge to make the method practical in geotechnical laboratories.

The table top centrifuge

The centrifuge used for the experiments is a Laboratory bench-scale centrifuge type Sigma 3–18, which costs less than 10.000 Euro. This set-up is equipped with two microprocessors for the independent control of the rotor recognition and the overspeed signal. The speed is continuously controlled by these microprocessors with an accuracy of 1 rpm to a maximum of 4200 rpm. Our soil sample containers can, however, only withstand 2500 rpm. The duration of centrifugation, acceleration and deceleration curves can be set. The centrifuge is equipped with four buckets. Uneven loading of oppositely located buckets may lead to imbalance, in that case the drive automatically switches off and an imbalance warning message is displayed. The maximum radius from the central rotor to the base of the buckets is 171 mm.

The soil sample is housed in a custom designed fully closed system where a constant pressure prevails, see Figure 1. This closed system is constituted by a Duran glass tube (LSB inner tube, inner diameter 21.55 mm $\times$ height 85 mm) containing the soil sample. This tube is screwed to a screw cap which is provided in the screw cap of an outlet reservoir (LSB centrifuge, glass thread GL 60, out diam. 62.3 mm $\times$ height 110 mm). The soil container tube is provided with a Borosilicate glass 3.3 Vitra POR filter welded to the tube with a thickness about 2.55 mm to 3 mm, a diameter of about 21.55 mm, and a nominal pore size 1 $\mathrm{\mu }$m–1.6 $\mathrm{\mu }$m (class 5, ISO/4793 = P1,6) for the tests done on kaolin clay and a nominal pore size 16 $\mathrm{\mu }$m–40 $\mathrm{\mu }$m (class 3, ISO/4793 = P40) for the tests done on the mixtures of kaolin and sand. The glass containers are provided of a sufficiently thick wall (3 mm) resistant to centrifugal forces.

Fig. 1 Centrifuge holder with sample.

Measurements

Sample and water chamber properties

The necessary length sizes (e.g. height of the soil) were measured using a Vernier caliper (accuracy 0.01 mm). The mass of the soil, inlet and outlet water were measured by means of laboratory scales (accuracy 0.01 g). Two types of samples are considered.

The first sample is a commercial processed kaolin Rotoclay HB (Goonvean, St. Austell, UK). This kaolin clay was chosen as reference material because it has been largely used in previous research.[20] Deionized water, produced using a water purification system, was used to prepare the samples and as reference permeant solution. The Electrical Conductivity of deionized water was EC $\le$ 3.9 $\mathrm{\mu }$S/cm and the pH was about 7.6. The samples are then prepared as a mixture of kaolin clay and deionized water. The amount of water in the mix was set to two times the liquid limit of the clay to obtain a slurry (w = 2 $\times$ 59$\%$ = 118$\%$). In order to create suitable reconstituted samples,([21]) the clayey soil and the water were mixed in a dough mixer for 15 min until a homogeneous distribution was observed. Next, the slurry was poured in the glass centrifuge tubes using a vibratory table to remove air bubbles. Then, the sample was housed in an outflow glass reservoir. These samples are then consolidated outside or inside the centrifuge.

The second sample is a mixture of kaolin and Mol sand. This was chosen as a more permeable test, and because comparison is possible with the investigation done in [22]. Mixtures of kaolin and clay have been prepared as follows: both clay and sand were dried in the oven at 105${}^{\circ }$ for 16 h. With a ceramic pestle, some clods of clay were crushed in a ceramic mortar. An amount of 90 g of dry Mol sand with 10 g of dry kaolin clay was mixed with a spoon thoroughly. The soil was then placed in the glass tubes and water was poured from the top to saturate the samples under 1 g conditions until the visible water front reached the base filter. The initial water content of the mixture was about 20$\%$.

The porosity of these soils can be deduced by considering the dry weight of the soil ($W_S$), total weight after obtaining saturated flow through the sample, ($W_T$), and the total volume of the consolidated samples, $V_T$. We obtain for the saturated porosity ${\mathit{\theta }}_s$,(1) $$\begin{array}{c}\hfill {\mathit{\theta }}_s=\frac{W_T-W_S}{{\mathit{\rho }}_w{V}_T},\end{array}$$(1) where ${\mathit{\rho }}_w$ is the density of water. This is the effective saturated porosity. Using the density of the soil, a larger porosity would be obtained; however, for the experiments, we must assume those extra voids to be not accessible (e.g. air bubbles locked inside soil particles).

Sample centre of gravity (COG) and water mass measurements

At the moment, the centrifuge must be stopped to perform measurements. Doing measurements takes around 10 min. The water in the outflow chamber is poured out and weighted, as is the water in the inflow chamber. This gives us the needed water mass measurements. For the determination of the COG of the soil in the tubes, two digital scale professional mini balances were used. The tube is placed horizontally on two sharp holders placed on the balances. The weight of the tube at these two points is recorded and the COG is calculated by means of the equilibrium of forces.

Mathematical model

The centrifuge

We consider a one-dimensional simplification of the sample. In this, the sample starts (left boundary) at the distance $r=r_0$ from the centre of the centrifuge and ends at the distance $r_L=r_0+L$, see Figure 2. A water chamber is placed at position $r\in (r_0-{\ell }_0,r_0)$ of the centrifuge. A second water chamber is placed at the outflow side, at $r\in (r_0+L,r_0+L+T)$. We indicate the height of the water in the inflow chamber as $\ell (t)$, and the height of the water in the outflow chamber as l(t). The rotational speed of the centrifuge given in radians per second is denoted by $\mathit{\omega }(t)$.

Fig. 2 Schematic of the centrifuge parts.

We never consider a completely dry sample, so a cross section of the sample can be in two states: partially saturated characterized by capillary flow and saturated characterized by pores which are completely filled by water.

The piesometric head h is defined from the pressure p as ${\mathit{\rho }}_w\mathit{gh}(r)=p(r)$. In a gravitational field, head is the height of the water column at point r giving rise to the given pressure. Within a centrifuge, this analogy is no longer present, but it is nevertheless custom to use h(r). In a centrifuge, a water level of $\ell (t)$ in the inflow chamber will give rise to a pressure at the top of the sample of(2) $$\begin{array}{c}\hfill p(r_0)={\int }_{r_0-\ell (t)}^{r_0}{\mathit{\rho }}_w{\mathit{\omega }}^{2}r\mathit{dr}=\frac{{\mathit{\rho }}_w{\mathit{\omega }}^2}{2}\ell (t)(2r_0-\ell (t)).\end{array}$$(2)  

Consolidation

As the samples are created in the laboratory, they will undergo consolidation when stresses are applied. We need to be careful to avoid stresses in the centrifuge that are higher than the stresses that occur naturally in the subsurface (we allow depth up to 50 m). On the other hand, the model under consideration will not contain consolidation, so the samples must be prepared and consolidated to the correct stresses before experiments are performed to determine the soil parameters.

The stress in soil mechanics that is important here is the effective stress, on the soil particles, given by ${\mathit{\sigma }}^{\prime }=\mathit{\sigma }-u$ (when u is positive), with $\mathit{\sigma }$ the total stress, and u the pore pressure. As we will do unsaturated experiments, we will obtain a negative pore pressure. When the pore pressure is nil, all weights must be carried by the soil particles, and ${\mathit{\sigma }}^{\prime }=\mathit{\sigma }$. Consider the worst case scenario, a sample with void spaces filled for 99.999% of water, which rotates at M rounds per minute. The total stress at the base of the sample in this case is in simplified form given by(3) $$\begin{array}{c}\hfill {\mathit{\sigma }}_b={\int }_{r_0}^{r_0+L}{\left(\frac{2\mathit{\pi }M}{60}\right)}^2({\mathit{\rho }}_w{\mathit{\theta }}_s+{\mathit{\rho }}_s(1-{\mathit{\theta }}_s))r\mathit{dr}={\left(\frac{2\mathit{\pi }M}{60}\right)}^2{\mathit{\rho }}_t({(r_0+L)}^2-r_0^2),\end{array}$$(3) where ${\mathit{\rho }}_t$ is the total density for the soil sample filled with water. We can compare this with stress at a depth D for the same type of soil which is saturated with water but with water table of 0 at this depth, using the formula ${\mathit{\sigma }}_D={\mathit{\rho }}_t\mathit{gD}$, so with ${\mathit{\sigma }}_b$ we can identify the depth(4) $$\begin{array}{c}\hfill D=\frac{{\mathit{\sigma }}_b}{{\mathit{\rho }}_{t}g}.\end{array}$$(4) In Table 1, typical results are given for the table top centrifuge under consideration, so $r_0=105$ mm, $L=40$ mm, ${\mathit{\rho }}_t=1.65$ g/cm${}^3$.

Table 1 Rotational speed and corresponding total stress at base for a sample with density ${\mathit{\rho }}_{t}1.65$ g/cm${}^3$, $r_0=105$ mm, $L=40$ mm. Corresponding depth in normal gravity, corresponding g-level at the midpoint of the sample (given by $N=R{\mathit{\omega }}^2/g$, with R the distance from the rotation center), and corresponding stress at that same midpoint are also given.

Speed	Total stress base ${\mathit{\sigma }}_b$	Depth D	g-level N at midpoint	Total stress midpoint
rpm	kPa	m	$-$	kPa
100	1.78	0.11	1.4	0.82
200	7.11	0.45	5.6	3.27
300	15.99	1.01	12.6	7.35
600	63.96	4.02	50.3	29.42
900	143.90	9.05	113.2	66.19
1200	255.82	16.10	201.2	117.68
1500	399.72	25.15	314.4	183.87
2400	1023.28	64.39	804.9	470.71

We see that speeds up to 2400 rpm are allowed, as they correspond in worst case assumption with a depth of 64 m. The stresses in the sample varied rapidly, indicating that consolidation done in the centrifuge can lead to inhomogeneous porosity over the sample. Note further that during normal consolidation, the sample is saturated, so the effective stress will be much lower than the given total stress. However, as soon as drainage starts, the effective stress will be close to the total stress from Table 1, and compression can occur in the sample. Experimentally, the biggest challenge is hence to consolidate the sample in such a way before the start of the experiments, that no further, or neglectable, compression of the sample occurs during the measurements. Fortunately, changes in height of the sample in the centrifuge are easily recognized, so it is simple to verify this requirement.

As comparison, the samples in [14] underwent centrifugation up to 1500 kPa to determine the soil retention curve from equilibrium curves, and up to 30% in height change between 80 kPa and 1500 kPa was observed, though the recovered soil retention curve matched the pressure plate result.

Saturated flow

Saturated flow ($h>0$) in a sample of fixed length L is determined by Darcy flow, which gives(5) $$\begin{array}{c}\hfill K_s{\mathit{\partial }}_r[{\mathit{\partial }}_{r}h-\frac{\mathit{\omega }{(t)}^2}{g}r]=0,r_0<r<r_0+L,\end{array}$$(5) where $K_s$ is the saturated hydraulic conductivity, with initial condition $h(r,0)=0,\mathit{\omega }(0)=0$ and boundary conditions at the top (left),(6) $$\begin{array}{c}\hfill h(r_0,t)=\frac{{\mathit{\omega }}^2}{2g}\ell (t)(2r_0-\ell (t)),\end{array}$$(6) while for the outflow at the bottom we can consider a water table height $H_B$(7) $$\begin{array}{c}\hfill h(r_0+L,t)=h_B=\frac{{\mathit{\omega }}^2}{2g}{H}_B(2(r_0+L)-H_B),\end{array}$$(7) which is provided by a filled outflow chamber with overflow lip (called the overflow cup), or free outflow, which results instead in the boundary condition(8) $$\begin{array}{c}\hfill h(r_0+L,t)=0.\end{array}$$(8) Equation (8(8) $$\begin{array}{c}\hfill h(r_0+L,t)=0.\end{array}$$(8) ) is due to the fact that excess water must be present at the outflow boundary before droplets can form at the bottom and detach.

One would expect that for the experiments, free outflow is (8(8) $$\begin{array}{c}\hfill h(r_0+L,t)=0.\end{array}$$(8) ) sufficient to obtain results. It is technically also easier than using an overflow cup. However, in our unsaturated experiments, problems are encountered when using free outflow. We theorize that this is because the used bottom filter allows air penetration, and during the deceleration faze, air penetrates the bottom part of the sample, greatly increasing the effective stress at the bottom of the sample, and leading to compression at the base and reduced flow. Using a water table at the bottom allows to ensure no air penetration is possible from the base, instead water from the overflow will be sucked into the sample, as we will show in Section 4.5. For future experiments, we plan to or use filters that allow no air penetration, or use (7(7) $$\begin{array}{c}\hfill h(r_0+L,t)=h_B=\frac{{\mathit{\omega }}^2}{2g}{H}_B(2(r_0+L)-H_B),\end{array}$$(7) ), so $H_B>0$, or use continuous measurements in the centrifuge which means no deceleration is needed to perform measurements and problems due to air entry at the bottom are avoided.

In practice, to obtain (7(7) $$\begin{array}{c}\hfill h(r_0+L,t)=h_B=\frac{{\mathit{\omega }}^2}{2g}{H}_B(2(r_0+L)-H_B),\end{array}$$(7) ) we use two outflow chambers, a smaller one with overflow to impose the boundary condition, and a larger one collecting the water that flows out of the overflow and is used in the inverse determination of the parameters.

Solving (5(5) $$\begin{array}{c}\hfill K_s{\mathit{\partial }}_r[{\mathit{\partial }}_{r}h-\frac{\mathit{\omega }{(t)}^2}{g}r]=0,r_0<r<r_0+L,\end{array}$$(5) ) under the given initial and boundary conditions for the unknown flux and $\ell (t)$ is straightforward.[23] With an inverse method, this allows to optimize the measurements of outflow water in terms of the saturated conductivity $K_s$, which is the only soil parameter present in the model. The result can be compared with the classical formula for hydraulic conductivity from two outflow measurements, see e.g. [24],(9) $$\begin{array}{c}\hfill K_s=\frac{L}{N(t_2-t_1)}\ln \frac{(H_1+L)(2r_0+L-H_2)}{(H_2+L)(2r_0+L-H_1)},\end{array}$$(9) where N is defined as $N=(r_0+L)\frac{{\mathit{\omega }}^2}{g}$, so it is the g-level at the base of the sample.

In the case of water table height $H_B$ at the base, we derived the analogue formula as(10) $$\begin{array}{c}\hfill K_s=\frac{L}{(r_0+L-H_B)\frac{{\mathit{\omega }}^2}{g}(t_2-t_1)}\ln \frac{(H_1+L-H_B)(2r_0+L-H_2-H_B)}{(H_2+L-H_B)(2r_0+L-H_1-H_B)}.\end{array}$$(10)  

Unsaturated flow

When considering unsaturated flow, Richards’ equation must be solved, given by(11) $$\begin{array}{c}\hfill \{\begin{array}{c}\hfill {\mathit{\partial }}_{r}q=0,q=-K_s[{\mathit{\partial }}_{r}h-\frac{\mathit{\omega }{(t)}^2}{g}r],h\ge 0,\\ \hfill {\mathit{\partial }}_t\mathit{\theta }=K_s{\mathit{\partial }}_r[k(\mathit{\theta })({\mathit{\partial }}_{r}h-\frac{\mathit{\omega }{(t)}^2}{g}r)],h<0,\end{array}\end{array}$$(11) where q is the flow flux, $\mathit{\theta }$ is the saturation of the porous media, $\mathit{\omega }(t)$ the angular speed of rotation (in radians per second), $K_s$ the saturated hydraulic conductivity, g the gravitational constant, function $K_{s}k(\mathit{\theta })$ the hydraulic conductivity in the unsaturated region and $k(\mathit{\theta })$ is called the relative permeability. The initial condition for an experimental run (after measurements) is the original saturation $\mathit{\theta }$ of the sample where unsaturated and $h(r,0)=0$ otherwise (as $\mathit{\omega }(0)=0$). The boundary condition at the top is (6(6) $$\begin{array}{c}\hfill h(r_0,t)=\frac{{\mathit{\omega }}^2}{2g}\ell (t)(2r_0-\ell (t)),\end{array}$$(6) ) as in saturated flow when water is present in the inflow chamber, and a no-flow boundary condition otherwise. The very first experiment is always saturated flow, hence unsaturated flow is always unsaturated at the top with possibly a saturated zone at the base.

At the base, we again need to distinguish between an overflow chamber or free outflow. In the last case, we have a no-flow boundary condition if $h(r_0+L)<0$, and otherwise $h(r_0+L)=0$. If, however, the bottom is kept at a known constant water level (overflow), the bottom will remain saturated, and we have $h(r_0+L)=h_B$ as in (7(7) $$\begin{array}{c}\hfill h(r_0+L,t)=h_B=\frac{{\mathit{\omega }}^2}{2g}{H}_B(2(r_0+L)-H_B),\end{array}$$(7) ), with a saturated zone on top of this boundary, separated from the unsaturated zone at a moving interface, which we denote by $r_s(t)$. As the samples are thin, we neglect the occurrence of a drying shock.[12]

Denote by $S_e=\frac{\mathit{\theta }-{\mathit{\theta }}_r}{{\mathit{\theta }}_s-{\mathit{\theta }}_r}$ the effective saturation where ${\mathit{\theta }}_s$ is the volumetric water content at saturation and ${\mathit{\theta }}_r$ the residual volumetric water content. We have $S_e\in (0,1)$, since $\mathit{\theta }\in ({\mathit{\theta }}_r,{\mathit{\theta }}_s)$. We consider soil hydraulic properties which were proposed by [25](12) $$\begin{array}{c}\hfill \begin{array}{cc}\hfill S_e& =\frac{1}{{(1+{\left(\mathit{\gamma }h\right)}^n)}^m},h\in (-\infty ,0),\hfill \\ \hfill k(S_e)& =S_e^{1/2}{[1-{(1-S_e^{1/m})}^m]}^2,\hfill \end{array}\end{array}$$(12) where $m=1-1/n$, $n>1$ and $\mathit{\gamma }$ are empirical soil parameters. Hence, flow in an unsaturated region can be rewritten in terms of hydraulic pressure head as(13) $$\begin{array}{c}\hfill {\mathit{\partial }}_h{S}_e{\mathit{\partial }}_{t}h={\mathit{\partial }}_r[K_{s}k(h)({\mathit{\partial }}_{r}h-\frac{{\mathit{\omega }}^2}{g}r)]\end{array}$$(13) where $k(h)$ is $k(S_e(h))$ from (12(12) $$\begin{array}{c}\hfill \begin{array}{cc}\hfill S_e& =\frac{1}{{(1+{\left(\mathit{\gamma }h\right)}^n)}^m},h\in (-\infty ,0),\hfill \\ \hfill k(S_e)& =S_e^{1/2}{[1-{(1-S_e^{1/m})}^m]}^2,\hfill \end{array}\end{array}$$(12) ).

Several numerical solution methods have been proposed to solve (13(13) $$\begin{array}{c}\hfill {\mathit{\partial }}_h{S}_e{\mathit{\partial }}_{t}h={\mathit{\partial }}_r[K_{s}k(h)({\mathit{\partial }}_{r}h-\frac{{\mathit{\omega }}^2}{g}r)]\end{array}$$(13) ). A common approach is based on using Picard iterations to overcome the non-linearities present in the mode.[26] This approach is applied in the often used HYDRUS code.[27] Another possibility is to leverage the available high accuracy ODE solvers via the method of lines, which was shown to work excellent for Richards’ equation in [28] using a mass-pressure approach. The same good result can be obtained via the method of lines using a saturation formulation or a mixed saturation-pressure formulation, augmented with a global mass balance equation,[19] which is the approach we also follow in this paper, due to its ease of implementation for 1D problems.

So, to solve (13(13) $$\begin{array}{c}\hfill {\mathit{\partial }}_h{S}_e{\mathit{\partial }}_{t}h={\mathit{\partial }}_r[K_{s}k(h)({\mathit{\partial }}_{r}h-\frac{{\mathit{\omega }}^2}{g}r)]\end{array}$$(13) ), the method of lines is used. In this, a finite volume discretization based on the flux $q$ is done in space on a discrete grid over the domain $r\in (r_0,r_0+L)$. This leads to a system of ODE, which is solved with the Sundials solver, see [23]. Note that the problem is singular where the transition from saturated to unsaturated occurs. To avoid problems, the interface between the saturated and unsaturated part is modelled, effectively splitting the domain in two parts that connect at the interface $r_s$, as was done in [19].

The soil parameters occurring in this model are now ${\mathit{\theta }}_s$, ${\mathit{\theta }}_r$, $K_s$ and the properties of the soil retention curve as proposed by van Genuchten, $\mathit{\gamma }$ and $n$. The effective porosity ${\mathit{\theta }}_s$ can be calculated from (1(1) $$\begin{array}{c}\hfill {\mathit{\theta }}_s=\frac{W_T-W_S}{{\mathit{\rho }}_w{V}_T},\end{array}$$(1) ), while $K_s$ is determined from the saturated experiments. The value of ${\mathit{\theta }}_r$ is only important in the cases it is possible to go to low saturation values. Before constructing the inverse method to determine the parameters of the soil retention curve, we present in details the solution method for the case of a water table at the base of the sample, as this has not been modelled before.

In details: drainage into an overflow cup

Drainage in an overflow leads to a partial differential equation with a dynamic boundary condition. At the start, and while the centrifuge rotates, water is drained from the sample, see Figure 3.

Fig. 3 Schematic representation of drainage into an overflow cup. Left: as the water drains, the height $H_B$ is constant. The interface s between partially saturated and saturated zone moves downward. Right: if the centrifuge slows down, part of the water from the overflow cup is sucked into the sample, and the interface s moves upward.

Denote $H_B^{\max }$ the maximal height of water in the overflow cup, $A_o$ the water surface area of the overflow cup, A the area of the sample and $r_B=r_L-H_B$ the distance from centrifuge axis to the basin water table. The model for drainage in an overflow cup is given by(14) $$\begin{array}{cc}\hfill {\mathit{\partial }}_t\mathit{\theta }& =K_s{\mathit{\partial }}_r[k(\mathit{\theta })({\mathit{\partial }}_{r}h-\frac{\mathit{\omega }{\left(t\right)}^2}{g}r)],r_0\le r<r_s,\hfill \end{array}$$(14) (15) $$\begin{array}{cc}\hfill q_{\mathit{sat}}& =-\frac{K_s}{L-s}\frac{\mathit{\omega }{\left(t\right)}^2}{2g}[r_s^2-r_B^2],\hfill \end{array}$$(15) (16) $$\begin{array}{cc}\hfill {\stackrel{·}{H}}_B(t)& =\{\begin{array}{cc}\frac{q_{\mathit{sat}}}{A_o}A\hfill & H_B<H_B^{\max }\hfill \\ \frac{q_{\mathit{sat}}}{A_o}A\hfill & H_B=H_B^{\max }\wedge q_{\mathit{sat}}<0\hfill \\ 0\hfill & H_B=H_B^{\max }\wedge q_{\mathit{sat}}\ge 0,\hfill \end{array}\hfill \end{array}$$(16) (17) $$\begin{array}{cc}\hfill {\stackrel{·}{M}}_{\mathrm{out}}\left(t\right)& =\{\begin{array}{cc}{\mathit{\rho }}_w{q}_{\mathit{sat}}A\hfill & q_{\mathit{sat}}\ge 0\wedge H_B=H_B^{\max }\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}\hfill \end{array}$$(17) with boundary conditions:(18) $$\begin{array}{cc}\hfill q\left(r_0\right)& =0,\hfill \end{array}$$(18) (19) $$\begin{array}{cc}\hfill h\left(r_L\right)& =\frac{{\mathit{\omega }}^2}{2g}(r_L^2-r_B^2),\hfill \end{array}$$(19) (20) $$\begin{array}{cc}\hfill h\left(r_s\right)& =0,\hfill \end{array}$$(20) (21) $$\begin{array}{cc}\hfill H_B\left(0\right)& =H_B^{\max },\hfill \end{array}$$(21) (22) $$\begin{array}{cc}\hfill M_{\mathrm{out}}\left(0\right)& =0.\hfill \end{array}$$(22) Equation (14(14) $$\begin{array}{cc}\hfill {\mathit{\partial }}_t\mathit{\theta }& =K_s{\mathit{\partial }}_r[k(\mathit{\theta })({\mathit{\partial }}_{r}h-\frac{\mathit{\omega }{\left(t\right)}^2}{g}r)],r_0\le r<r_s,\hfill \end{array}$$(14) ) describes the saturation evolution in the unsaturated part, the saturated flow formula (15(15) $$\begin{array}{cc}\hfill q_{\mathit{sat}}& =-\frac{K_s}{L-s}\frac{\mathit{\omega }{\left(t\right)}^2}{2g}[r_s^2-r_B^2],\hfill \end{array}$$(15) ) is derived directly from Darcy’s flow in the saturated part (at the bottom) and determines the water height table $H_B(t)$ in the overflow chamber described by (16(16) $$\begin{array}{cc}\hfill {\stackrel{·}{H}}_B(t)& =\{\begin{array}{cc}\frac{q_{\mathit{sat}}}{A_o}A\hfill & H_B<H_B^{\max }\hfill \\ \frac{q_{\mathit{sat}}}{A_o}A\hfill & H_B=H_B^{\max }\wedge q_{\mathit{sat}}<0\hfill \\ 0\hfill & H_B=H_B^{\max }\wedge q_{\mathit{sat}}\ge 0,\hfill \end{array}\hfill \end{array}$$(16) ), (21(21) $$\begin{array}{cc}\hfill H_B\left(0\right)& =H_B^{\max },\hfill \end{array}$$(21) ) and the amount of expelled water described by (17(17) $$\begin{array}{cc}\hfill {\stackrel{·}{M}}_{\mathrm{out}}\left(t\right)& =\{\begin{array}{cc}{\mathit{\rho }}_w{q}_{\mathit{sat}}A\hfill & q_{\mathit{sat}}\ge 0\wedge H_B=H_B^{\max }\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}\hfill \end{array}$$(17) ), (22(22) $$\begin{array}{cc}\hfill M_{\mathrm{out}}\left(0\right)& =0.\hfill \end{array}$$(22) ). Boundary condition (18(18) $$\begin{array}{cc}\hfill q\left(r_0\right)& =0,\hfill \end{array}$$(18) ) is the no-flow condition at $r_0$ and can be equally written as ${[{\mathit{\partial }}_{r}h-\frac{{\mathit{\omega }}^2}{g}r]}_{r=r_0}=0$. Condition (20(20) $$\begin{array}{cc}\hfill h\left(r_s\right)& =0,\hfill \end{array}$$(20) ) indirectly determines the position of the interface $r_s$, which divides the saturated zone from the unsaturated. Finally, (19(19) $$\begin{array}{cc}\hfill h\left(r_L\right)& =\frac{{\mathit{\omega }}^2}{2g}(r_L^2-r_B^2),\hfill \end{array}$$(19) ) is the boundary condition for hydraulic pressure at the sample end $r_L$, which is implied by the hydrostatic pressure of water in the overflow basin. During normal operation, we have $H_B(t)=H_B(0)$ and the interface moves towards $r_0+L-H_B$. However, when the centrifuge slows down, the pressure at the base of the sample drops, and water is sucked into the sample. This causes $H_B(t)$ to drop, and the interface $r_s$ to move towards $r_0$.

During operation, the system moves towards an equilibrium. For this setting, the interface will move faster to its equilibrium position than the saturation profile. Numerically, it is important to avoid instability due to the interface being close to its equilibrium position. This position follows from following lemma.

Lemma 4.1

If during sample drainage (constant $\mathit{\omega }$, initial $q_{\mathit{sat}}(t_0)>0$) it holds $r_s\left(t_0\right)\le r_B\left(t_0\right)$, then this holds also $\forall t\ge t_0$.

 

Proof

The pressure at the end of the soil sample is caused by the hydrostatic pressure of water in the outflow chamber. From Darcy’s equation, it follows that the flux through the saturated part of the sample is given by(23) $$\begin{array}{c}\hfill q_{\mathit{sat}}=-\frac{K_s}{L-s}\frac{{\mathit{\omega }}^2}{2g}[r_s^2-r_B^2],\end{array}$$(23) and hence,(24) $$\begin{array}{c}\hfill \frac{d{q}_{\mathit{sat}}}{\mathit{dt}}=-\frac{K_s}{L-s}\frac{2\mathit{\omega }}{2g}\frac{d\mathit{\omega }}{\mathit{dt}}[r_s^2-r_B^2]-\frac{K_s}{L-s}\frac{{\mathit{\omega }}^2}{2g}2r_s\frac{{\mathit{dr}}_s}{\mathit{dt}}.\end{array}$$(24) During drainage, we have $\frac{d\mathit{\omega }}{\mathit{dt}}=0$, and from (23(23) $$\begin{array}{c}\hfill q_{\mathit{sat}}=-\frac{K_s}{L-s}\frac{{\mathit{\omega }}^2}{2g}[r_s^2-r_B^2],\end{array}$$(23) ) it follows that when the saturated front reaches the water level in the outflow chamber, the saturated flux becomes zero. Therefore, the saturated front can from then onward only change due to flux from the unsaturated part, and will hence not continue moving forward. So $\frac{{\mathit{dr}}_s}{\mathit{dt}}\le 0$, $q_{\mathit{sat}}\ge 0$. Hence, $r_s=r_0+s\le r_0+L-H_B=r_B$. $\square$

As a consequence of this Lemma, if the interface is within the numerical precision from $r_B$, we can recast the problem to one with a fixed interface at $r_B$. To determine the position of the moving interface, we use an algebraic equation based on total mass balance of the water held in the entire experiment (like in [19]). So, this includes the water contained in the sample, water inside the basin and expelled water captured in the outflow chamber. The resulting system is then a system of differential algebraic equations (DAE), which requires good initial guesses for the time derivatives. Typically, a DAE solver may have problems finding the initial value of $\stackrel{·}{s}$, in this case, providing as initial guess $\stackrel{·}{s}=0.1$ was sufficient.

Important in the calculation of a solution with a DAE solver is the creation of a suitable space grid. We use a discretization so that the points $r_i$ create a geometric sequence, i.e. if $\mathrm{\Delta }{r}_i=(r_{i+1}-r_i)$ then $\mathrm{\Delta }{r}_i=k\mathrm{\Delta }{r}_{i-1}=k^i\mathrm{\Delta }{r}_0,i>0$. Therefore, $L=(1+k+k^2+\dots +k^N)r_0$, where $N$ is the number of discretization points which together with the coefficient $k$ is given as user input, and allows to compute $r_0$. We typically use $k\in 〈0.90,0.94〉$, which refines the grid towards the saturated interface $r_s$.

Inverse determination soil retention curve

The numerical solution $S_e$ of the unsaturated flow results in a computed value of the effective saturation in terms of the applied parameters $p=(\mathit{\gamma },n)$, in other words, we obtained $S_e(r,t|\mathit{\gamma },n)$. We also obtain the inflow flux, and the outflow flux, by keeping track of the boundary conditions. Based on the computed $S_e$, we can determine the COG of the water, denoted by $G_w(t,\mathit{\gamma },n)$, and from the boundary conditions, the amount of inflow water and outflow water, denoted by $M_{\mathit{wi}}(t,\mathit{\gamma },n)$ and $M_{\mathit{wo}}(t,\mathit{\gamma },n)$, respectively. Note that all these measurements are global data; that is, they are independent on the position r.

The computed values should be as close as possible to the measured values at discrete times $t_1,\dots ,t_m$. We indicate the measured values with a bar, e.g. ${\overline{G}}_w(t_1)$. The cost functional to minimize in terms of the parameters is hence(25) $$\begin{array}{cc}\hfill \mathcal{F}(\mathit{\gamma },n)& =w_1\sum _{i=1}^m{(G_w(t_i,\mathit{\gamma },n)-{\overline{G}}_w(t_i))}^2+w_2\sum _{i=1}^m{(M_{\mathit{wi}}(t,\mathit{\gamma },n)-{\overline{M}}_{\mathit{wi}}(t_i))}^2\hfill \\ \hfill & +w_3\sum _{i=1}^m{(M_{\mathit{wo}}(t,\mathit{\gamma },n)-{\overline{M}}_{\mathit{wo}}(t_i))}^2,\hfill \end{array}$$(25) where $w_1,w_2,w_3$ are weights given to the different types of measurements. For the minimization of this cost functional, the Levenberg–Marquardt method (LMM) is used, which gives a solution $({\mathit{\gamma }}^{\ast },n^{\ast })$. If the value of the cost functional for this solution is sufficiently low, we assume a global minimum is found. Many free implementations of LMM are available, we use the lmdif routine of the minpack library. However, several technical difficulties arise for the considered problem, that are handled in the following way:

1. Many possible parameters lead to extremely slow initial condition calculations of the DAE system. To alleviate this, we use several tricks: an initial interface $r_s$ is introduced at $r_0+0.8$ cm. When computing over a relatively wide range of van Genuchten parameters $\mathit{\gamma },n$, we use a dynamic range of initial head $h_{\mathit{init}}$ in the unsaturated zone. We choose $h_{\mathit{init}}$ such that it is within a specified interval, which in our case is set to $〈-0.05,-5.0〉$. The reason is that if $\left(\mathit{\gamma }h\right)$ appearing in the van Genuchten formula $S_e=\frac{1}{{(1+{\left(\mathit{\gamma }h\right)}^n)}^m}$ is small, the relative saturation $u$ is close to $1$, which can cause a large slow down of the computation of derivatives and interface $r_s$. Therefore, we specify a constant value $c_{\mathit{init}}$ such that $\left(\mathit{\gamma }{h}_{\mathit{init}}\right)=c_{\mathit{init}}$ from which we obtain the $h_{\mathit{init}}$ value. We use this value if it falls into our specified range, otherwise the closest edge value of the interval is taken. Knowing the estimates of the optimal values, this dynamic setting of $h_{\mathit{init}}$ can be turned off so as to investigate if it influences the final result. Finally, at the interface we set $h(s)=-0.05$ instead of the expected $h(s)=0$.

2. We specify ranges where parameters are to be expected. Also, as $K_s>0$ and/or $\mathit{\gamma }<0$, we use logarithmic transformation for these together with the user specified ranges.

3. Many possible $n$ and $\mathit{\gamma }$ lead to slow or impossible computations of the direct problem—or simply they fall outside of the expected range. To allow the LMM method to recover, we introduce fake costs for these parameters. In practice, those outside the interval are penalized based on the distance from the interval. In case the solver could not finish a computation inside the interval, we penalize with both distances to the edges. The penalization function used is $f_{\mathrm{pen}}(p)=e^{p-p_{\mathrm{boundary}}}+e^{-(p-p_{\mathrm{boundary}})}$ where $p_{\mathrm{boundary}}$ is the interval edge value used for the parameter.

In above, the unknown parameters are $p=(\mathit{\gamma },n)$, but this can be readily extended to $p=(\mathit{\gamma },n,K_s,{\mathit{\theta }}_r,{\mathit{\theta }}_s)$. Trying to recover five parameters from the available measurements will, however, not lead to good results. The best approach is to determines $K_s$ from a saturated test, and ${\mathit{\theta }}_s$ by oven-drying the sample.

Experiments

Saturated conductivity

The saturated conductivity was determined for the different filters, and the different soil samples prepared. The results are given in Figure 4 for the first four filters. Here, the results are computed with the model (5(5) $$\begin{array}{c}\hfill K_s{\mathit{\partial }}_r[{\mathit{\partial }}_{r}h-\frac{\mathit{\omega }{(t)}^2}{g}r]=0,r_0<r<r_0+L,\end{array}$$(5) ) where $K_S$ is optimized to match the measurements, which gives a difference of less than 1.5% compared with the classical formula (9(9) $$\begin{array}{c}\hfill K_s=\frac{L}{N(t_2-t_1)}\ln \frac{(H_1+L)(2r_0+L-H_2)}{(H_2+L)(2r_0+L-H_1)},\end{array}$$(9) ), so these lasts have not been added to the figure.

Fig. 4 Saturated permeability for the first 4 filters used.

Important to note here is that the measurements of inflow/outflow water in the centrifuge is sufficiently accurate to obtain good results that are reproducible. The variation in the centrifuge tests can be attributed to measurement uncertainty (duration rotation, correct rpm, startup/slowdown curve), while the difference in obtained permeability for the falling head test can be attributed to the difficulty to keep the filter saturated in the experimental set-up for the falling head test. The given falling head value is the average of two falling head tests, which themselves differ sometimes more than with the centrifuge result. A single falling head test for these filters took more than 150 min, while at 1000 rpm only 4 min of centrifugation time is required.

It is also important to note that in the centrifuge, free outflow is used, so no head value is imposed on the base of the filter. The centrifugal force is sufficiently high to prevent air penetration. The other two filters used were only tested with a normal falling head test. We find for filter 16 $K_S=1.351{10}^{-5}$ m/s, and for filter 17, $K_S=1.103{10}^{-5}$ m/s. These filters are used for more permeable samples.

This procedure has also been performed on the different soil samples with good results. In this case, the filter permeability is taken into account in the modelling, and the soil sample permeability is obtained. Considering a preconsolidated kaolin sample, we present in Table 2, next to the already mentioned results, the results of the classical falling head test in a normal gravity field, compared with the average of least squares optimizations (like (25(25) $$\begin{array}{cc}\hfill \mathcal{F}(\mathit{\gamma },n)& =w_1\sum _{i=1}^m{(G_w(t_i,\mathit{\gamma },n)-{\overline{G}}_w(t_i))}^2+w_2\sum _{i=1}^m{(M_{\mathit{wi}}(t,\mathit{\gamma },n)-{\overline{M}}_{\mathit{wi}}(t_i))}^2\hfill \\ \hfill & +w_3\sum _{i=1}^m{(M_{\mathit{wo}}(t,\mathit{\gamma },n)-{\overline{M}}_{\mathit{wo}}(t_i))}^2,\hfill \end{array}$$(25) ), minimize the difference between 10 centrifuge measurements of the outflow water and the model results) at different rpm, to determine $K_S$. Note that by using only outflow water, we can also include experiments that show some consolidation, as this is mainly visible in the volume of water of the inflow water chamber. Due to the pre-consolidation, consolidation will, however, be very limited.

Table 2 Permeability of the soil samples and filters.

Name	Falling head	Centrifuge model
	$K_S$ [m/s]	avg. $K_S$ [m/s]
Filter tube 1	$5.53{10}^{-8}$	$9.32{10}^{-8}$
Filter tube 2	$4.95{10}^{-8}$	$7.50{10}^{-8}$
Filter tube 4	$9.89{10}^{-8}$	$7.44{10}^{-8}$
Filter tube 5	$8.17{10}^{-8}$	$6.24{10}^{-8}$
Kaolin tube 1	$1.136{10}^{-8}$	$1.173{10}^{-8}$
Kaolin tube 2	$1.209{10}^{-8}$	$1.139{10}^{-8}$
Kaolin tube 4	$1.331{10}^{-8}$	$1.137{10}^{-8}$
Kaolin tube 5	$1.259{10}^{-8}$	$1.117{10}^{-8}$
Filter tube 16	$1.351{10}^{-5}$	–
Filter tube 17	$1.103{10}^{-5}$	–
Mixture tube 16	–	$2.096{10}^{-6}$

The last line in Table 2 shows the saturated conductivity for the kaolin-sand mixture based on centrifuge saturated flow. No falling head test was performed on this sample.

Unsaturated conductivity kaolin

For kaolin, the current centrifuge is not able to drain the water as this experiment could only attain 1000 rpm before shaking stopped the centrifuge automatically. Consolidation occurs with only some 3.6 mm of outflow water after 15 h, of which 2.4 can be attributed to drainage and 1.2 due to consolidation. This is too little information to use the model to determine the soil retention curve. The model should also be extended to contain consolidation to process this data.

The model can, however, be used to predict the flow behaviour at different rotational speeds. A pressure plate test was done to determine the soil parameters of the kaolin used, from which we obtained $\mathit{\gamma }=-0.0057$ cm${}^{-1}$, $n=1.36$, ${\mathit{\theta }}_s=0.6$, see Figure 5. This is combined with saturated conductivity $K_S=1.14{10}^{-8}$ m/s, with the results presented in Figure 6. There the evolution of the head and saturation can be seen for rotational speed 1000 rpm (as used), and 6000 rpm, over 60 min of operation.

Fig. 5 Measured soil retention curve, and interpolated value for the kaolin sample.

Fig. 6 Drainage of a kaolin sample in the centrifuge at $\mathit{\omega }=1000$ rpm (top) and $\mathit{\omega }=6000$ rpm (bottom) over 60 min starting from saturated.

We can conclude that a much higher rotational speed is needed to determine the full soil retention curve, but at 6000 rpm it would be possible to recover most of the curve. However, this might lead to consolidation pressures that are not present in the field. As mentioned, the currently used centrifuge starts to vibrate excessively, which does not allow us to try higher speeds and verify this. The model can, however, be used to simulate the expected behaviour. Note that the model has no problem with the fact that n is close to 1, a typical issue noted in [29]. As drainage with kaolin was not possible, we could not investigate if the adaptations suggested in [29] to the van Genuchten model are required. Instead, we focus next on drainage of sand–clay mixtures.

Unsaturated conductivity sand–clay mixture

The developed centrifuge does perform good for other type of soils which are more permeable than clays. To illustrate this, a drainage experiment was performed on the 100 kPa dry pre-loaded 90% sand–10% kaolin mixture, which was used for the saturated conductivity test. The sample has a void ratio $e=0.538$, which corresponds to a porosity of ${\mathit{\theta }}_s=0.3498$. Typical values measured for expelled water and gravitational centre (GC) are given in Table 3. Here, every line is a new measurement, after the given time in seconds. Performing a measurement takes on average 10 min outside of the centrifuge. The entire measurement was done in one day (validating the measurement techniques and procedures, and constructing the mathematical model took up most of the research time, but these are one time setup costs, which need not be repeated).

Table 3 Drainage experiment on 90% sand-10% kaolin mixture. Every line is a new measurement, after the given time in seconds.

Time (s)	rpm	Height soil	Height expelled	Cumulative	Grav. Center
		at start (cm)	water (cm)	height (cm)	(cm)
60	600	5.731	0.03896	0.03896	2.525
120	600	5.701	0.04174	0.08070	2.483
240	600	5.701	0.09461	0.17532	2.387
480	600	5.701	0.06679	0.24211	2.255
960	600	5.701	0.03896	0.28107	2.154
1920	600	5.701	0.03061	0.31168	2.085
3840	600	5.701	0.00278	0.31446	2.060
60	1200	5.701	0.07514	0.38960	1.997
120	1200	5.676	0.05566	0.44525	1.947
240	1200	5.676	0.05566	0.50091	1.886
480	1200	5.676	0.03896	0.53987	1.794
960	1200	5.676	0.05009	0.58996	1.686
1920	1200	5.676	0.01113	0.60109	1.671
60	2400	5.676	0.00557	0.60666	2.412
120	2400	5.456	0.01391	0.62057	1.639
240	2400	5.646	0.02505	0.64562	1.592
480	2400	5.646	0.02226	0.66788	1.485
960	2400	5.646	0.00835	0.67623	1.416
1920	2400	5.646	0	0.67623	1.392

The measurements are used in the inverse problem to determine the saturated conductivity $K_s$ and the van Genuchten parameters $\mathit{\gamma }$ and n. The overflow in this experiment is such that the water table was maintained around 2 cm above the base of the sample. We obtain $K_S=1.16{10}^{-7}$ m/s, $\mathit{\gamma }=-0.0136$ cm${}^{-1}$ and $n=1.96$ if we give the expelled water measurements double the weight of the GC ($w_1=1$, $w_2=0$, $w_3=2$ in (25(25) $$\begin{array}{cc}\hfill \mathcal{F}(\mathit{\gamma },n)& =w_1\sum _{i=1}^m{(G_w(t_i,\mathit{\gamma },n)-{\overline{G}}_w(t_i))}^2+w_2\sum _{i=1}^m{(M_{\mathit{wi}}(t,\mathit{\gamma },n)-{\overline{M}}_{\mathit{wi}}(t_i))}^2\hfill \\ \hfill & +w_3\sum _{i=1}^m{(M_{\mathit{wo}}(t,\mathit{\gamma },n)-{\overline{M}}_{\mathit{wo}}(t_i))}^2,\hfill \end{array}$$(25) )). As we have GC and outflow water measurements, we can also obtain the values taking into account only specific measurements. If we only use GC, we obtain $K_S=3.1{10}^{-7}$ m/s, $\mathit{\gamma }=-0.0154$ cm${}^{-1}$ and $n=2.24$, while using only outflow water mass we obtain $K_S=3.8{10}^{-7}$ m/s, $\mathit{\gamma }=-0.0234$ cm${}^{-1}$ and $n=1.44$. The experiments and their best fit by the three cases is given in Figure 7. The GC and outflow results there can be seen as the extremal curves, with depending on the weight given to the experiments, the combined result lying in between. As can be expected, GC matches closest the GC measurements, but has largest error for outflow. This validates the use of both measurement types.

Fig. 7 Experiments performed with the centrifuge. Crosses indicate the measured values, Triangles up the model results for the optimal value taking Gravitational center (GC) and Outflow (MO) into account. Triangles down only using GC, Triangles left using only MO. Full lines indicate the results using the pressure plate fitted soil parameters with two different values of $K_S$.

Fig. 8 Measured and inversely determined soil retention curves. Reference curve (Ref.) is the measurements done with ASTM D2325, with the other lines best measurement fit from the model.

The $\mathit{\gamma }$ and n values are compared with those obtained from a pressure plate measurement, which was performed by an outside laboratory according to ASTM D2325 on eight different samples. It took more than two months before the results were obtained. Although the actual experimental time will have been much less, it will have been considerably more than the couple of days needed with the bench scale centrifuge. Of the obtained result, one measurement was identified as erroneous, probably because too little time was given to reach equilibrium (a saturation value of 0.081 at 97885 Pa). Doing a least-squares fit with the van Genuchten model for the soil retention curve, we obtain as best fit $\mathit{\gamma }=-0.018\pm 0.006$ cm${}^{-1}$ and $n=2.18\pm 0.53$, which we combine once with $K_S=2.096{10}^{-6}$ m/s and once with $K_S=2.096{10}^{-7}$ m/s to compare with experiments in Figure 7 (labelled Ref., these give the theoretical discrepancy with the measurements as predicted by the mathematical model). The resulting soil retention curves are given in Figure 8. Observing this result, it is clear that using GC or a combination of GC and outflow gives a very good match with the result from the ASTM D2325 experiment. Using only MO gives a bad match, which is evident in the worse recovery of the n parameter. We also note that the recovered $K_S$ is an order less than the one determined on another sample of the same material in Table 2. As the numerical analysis was finished some months after the experiment, this could not be further investigated. Fixing $K_S$ in the inverse determination to the value obtained in the saturated test, $2.096{10}^{-6}$ does not give good results, as is evident in Figures 7 and 8, where results of the model with that $K_S$ have been added. We conclude that the sample used in the draining indeed had a lower saturated conductivity, which is consistent with the expelled water found. Independent of differences possible in saturated conductivity, the soil retention curve found matches that of the samples used in the ASTM D2325 experiment. An overview of all parameters, with corresponding determined deviation, is given in Table4.

Table 4 Resulting parameters arising in the model. If deviation is given, the parameter was determined via non-linear least-squares, if not given, it was known input.

Description	$K_s$ (m/s)	$\mathit{\gamma }$ (cm${}^{-1}$)	n	${\mathit{\theta }}_s$
Pressure plates (ASTM D2325)	–	$-0.0180\pm 0.006$	$2.18\pm 0.53$	$0.368\pm 0.035$
Based on outflow	$(3.8\pm 1.3){10}^{-7}$	$-0.0234\pm 0.0043$	$1.44\pm 0.05$	$0.346$
Based on GC	$(3.1\pm 0.5){10}^{-7}$	$-0.0154\pm 0.0011$	$2.24\pm 0.07$	$0.346$
Based on outflow & GC	$(2.0\pm 1.5){10}^{-7}$	$-0.0121\pm 0.0022$	$2.17\pm 0.37$	$0.346$
Based on outflow & 2 $\times$ GC	$(2.6\pm 0.9){10}^{-7}$	$-0.0136\pm 0.0014$	$2.26\pm 0.19$	$0.346$
Based on 2 $\times$ outflow & GC	$(1.16\pm 1.3){10}^{-7}$	$-0.0136\pm 0.0044$	$1.96\pm 0.40$	$0.346$
Outflow & GC, fixed $K_s$	$2.096{10}^{-6}$	$-0.0109\pm 0.0023$	$1.69\pm 0.07$	$0.346$
Outflow & GC, fixed $K_s$	$2.096{10}^{-7}$	$-0.0122\pm 0.0022$	$2.15\pm 0.13$	0.346
Outflow & GC, 1200 rpm	$(2.3\pm 0.9){10}^{-7}$	$-0.0107\pm 0.0014$	$2.39\pm 0.29$	0.346

The time evolution of the saturation in the centrifuge, as calculated by the model, is given in Figure 9. The saturated part at the base of the sample (right) is clearly visible. We see that the lowest saturation is around 0.3, which corresponds to a minimal water content of 0.11. Higher rotational speeds would be required to drain the sample more, and confidently determine the soil retention curve for lower water content. The current centrifuge 2400 rpm corresponds to a maximum capillary pressure of 30,000 Pa (compare with Figure 8). In Table 4, we add the recovered parameters should only the data up to 1200 rpm be used. The recovered parameters are immediately less good, though still usable. Specifically, the value of $\mathit{\gamma }$ is off, resulting in a lower estimated permeability than is the case in reality. We conclude that for an accurate recovery of the soil retention curve, it is needed to use an rpm that drains 80% of the sample in the unsaturated zone.

Fig. 9 Saturation curves as obtained by the model during operation of the centrifuge.

Overall, the experimental values can reasonably be recovered. The two parameters present in this model make for a sufficiently simplified model to obtain a good approximation from the experiments performed. From a geotechnical point of view, the centrifuge determined soil retention curve from outflow combined with GC can be used for predictive models.

The values can also be compared with the results in [22], which for such a mixture gives $n=1.20$, $\mathit{\gamma }=-0.028$ cm${}^{-1}$ to $-0.08$ cm${}^{-1}$, $K_s=1.5{10}^{-7}$ m/s for porosities (0.28–0.34). The $\mathit{\gamma }$ value differs a lot, indicating that our sample has a much higher bubbling pressure. The saturated conductivity from the unsaturated test does correspond with this value however. We can conclude that our sample is different from the one in [22] which is apparent also in the different void ratio. We attribute this difference to the different material preparation method. This does proof the danger of using literature data to set soil parameters. Doing an experiment to determine the soil parameters is very important, and the centrifuge method can provides these fast.

Conclusion

We conclude that based on transient data of outflow water and gravitational centre, the soil retention curve can be reconstructed, with a result matching the retention curve as found with an ASTM D2325 experiment.

We expect even better results to be obtained when measurements can be done while the centrifuge is rotating, which will be the next focus in this research. The fact that using only GC gave also good results is interesting. More tests should be performed to investigate if this is a general characteristic.

As the ASTM D2325 measurements allow for the possibility of different regimes in the sample, using a free-form determination for the soil retention curve instead of the van Genuchten model could also be considered. It remains to be seen if the amount of data is then sufficient for a reliable reconstruction. The code to analyse centrifuge experiments can be found at https://gitorious.org/centrifuge-1d/.

Acknowledgments

This research project is funded by IWT (Agency for Innovation by Science and Technology) and Geosound.be, within the project CENPERON (Centrifuge for the determination of the permeability of unsaturated soils), IWT project 110317.