The distribution of saturated clusters in wetted granular materials

The hydro-mechanical behaviour of partially saturated granular materials is greatly influenced by the spatial and temporal distribution of liquid within the media. The aim of this paper is to characterise the distribution of saturated clusters in granular materials using an optical imaging method under different water drainage conditions. A saturated cluster is formed when a liquid phase fully occupies the pore space between solid grains in a localized region. The samples considered here were prepared by vibrating mono-sized glass beads to form closely packed assemblies in a rectangular container. A range of drainage conditions were applied to the specimen by tilting the container and employing different flow rates, and the liquid pressure was recorded at different positions in the experimental cell. The formation of saturated clusters during the liquid withdrawal processes is governed by three competing mechanisms arising from viscous, capillary, and gravitational forces. When the flow rate is sufficiently large and the gravity component is sufficiently small, the viscous force tends to destabilize the liquid front leading to the formation of narrow fingers of saturated material. As the water channels along these liquid fingers break, saturated clusters are formed inside the specimen. Subsequently, a spatial and temporal distribution of saturated clusters can be observed. We investigated the resulting saturated cluster distribution as a function of flow rate and gravity to achieve a fundamental understanding of the formation and evolution of such clusters in partially saturated granular materials. This study serves as a bridge between pore-scale behavior and the overall hydro-mechanical characteristics in partially saturated soils.


Introduction
In soil mechanics, multiphase fluid flows in porous medium have been subjected to a considerable number of studies over the past decades since it is important to a broad range of applications, including oil recovery, risk assessment of landslides, and groundwater hydrology [1,2]. Strong variations exist during mixing two or more immiscible fluids in porous medium, depending on two dimensionless numbers, i.e., the Capillary number (Ca) and Bond number (Bo) [3,4]. Experimental and numerical studies concerning drainage processes in porous media often tend to focus on the evolution of viscous fingering at various flow rates [5,6]. However, studies on the morphological characteristics of liquid phase entrapment behind the front (interface between wetting and non-wetting fluids) remain relatively scarce.
This study aims to find a link between the two aforementioned dimensionless numbers and saturated cluster size distribution by establishing an experimental system where capillary, gravity and viscous forces can be examined. Here in, we focused on saturated clusters, which were formed when a liquid phase fully occupied the pore space between solid grains in a localised region. We have studied the saturated clusters obtained by draining water from a porous medium. To tune the gravitational and viscous forces, different tilting angles and drainage rates were employed on the system. The temporal and spatial information of the porous medium was characterised by means of a digital camera. The size distribution of saturated clusters were extracted from a series of images and statistically analysed by employing the maximum likelihood estimation [7].
The present paper is organised as follows. The experimental setup and methods for data analyses are described in Section 2. In Section 3, the saturated cluster distributions under various loading scenarios are presented and discussed. Finally, conclusions are presented in Section 4.

Experimental setup and methods
The experimental system was designed for temporally and spatially resolved imaging of interface motions in micro-models consisting of glass beads. The key elements of the experimental setup are shown in Fig.  1(a). A two-dimensional porous medium was constructed by a single layer of glass beads with a diameter of 2 mm between two parallel glass plates. The glass plates and beads were then slightly compressed by a few clamps, and silicone glue were applied on the bottom and side edges, leaving the top side connected to the atmosphere. The observable area was 180 mm 270 mm, and the resulting porosity was about ݊ ≈0.4. After the model was installed, the outlet channel was connected to a dispensing pump (Masterflex 77200-62, USA) that controlled a stable volumetric injection and withdrawal flow rate. A digital camera (Nikon D3300) controlled by a computer over digiCamControl software was used to record the experiments at given time eclipses based on different flow rates. Each image contains 1000 1500 pixels, which corresponds to a spatial resolution of 30 pixels per mm 2 . An electroluminescent panel provided sufficient background illumination for different shutter speeds. Both the light source and model could be tilted spontaneously to various inclination angles to vary the gravitational component on the fluid flows in the porous medium. The component of gravity along the model ݃' is given by ݃ sin(ߙ), where ݃ is the gravitational acceleration and ߙ is the inclination angle.
Water was used as the wetting fluid, dyed with dye tracer (Cole-Parmer 00298-06), and air was the nonwetting fluid invading from the open top edge. Experimental values of the wetting fluid are calculated by assuming room temperature at 25°C. The surface tension between two fluids ߛ is 72.75 mN/m. The wetting fluid has a viscosity ߤ of 1.00 mPa.s and a density ߩ of 998 kg/m3. The permeability of the medium ݇ was estimated from the Kozeny-Carman relation as 4.89 10-9 m 2 [8]. The applied flow rates ranged from 1 ml/min to 20 ml/min for different tilt angles. These resulted in the mean front velocity determined from volumetric flow rates at the cross-section area taking into account the porosity (shown in Table 1).
To quantify magnitudes of viscous or gravitational forces relative to capillary forces in the following experiments and analyses, we employed dimensionless Capillary (Ca) and Bond (Bo) numbers which are defined as follows [9]: with the wetting liquid density (ߩ), the component of gravitational acceleration ( ݃' ), typical pore size ( ܽ ) assumed to be 1/3 of bead diameter (more details can be found in [10]), viscosity of the wetting fluid (ߤ), surface tension (ߛ), the mean front velocity ‫,)ݒ(‬ permeability of the medium (݇).
In previous studies on the liquid front, it is found that the transition from stable to unstable front morphology conditions occurs at Bo*=0. The unstable front morphology exhibits capillary or viscous fingering during drainage stages. If the generalised Bond number is larger than zero, the gravitational force will overcome the viscous effects and lead to a flat plateau of the liquid front. However, when it is less than zero, the geometry of the front converts into unstable patterns.
The grayscale pixel-intensity histogram of raw images have two distinguishable peaks corresponding to the non-wetting and wetting fluids. In order to separate the air-filled region from the saturated region, all grey  levels under a given threshold were set to black ('wet') while the others were white ('dry'). The thresholds varied from images to images to account for the frequency changes of the electroluminescent panel. After the thresholding process, the whole observation area was separated into different regions with glass beads determined to be part of wet or dry region. This treatment did not affect the following analyses since only the saturated regions that are larger than the projection area of a typical glass bead (ߨ mm 2 in this study) were considered. After removal of trapped liquid volume with the projection area smaller than a single glass bead, the resulting black and white picture could provide the required geometrical information. To allow visualising separate clusters of different sizes, analyses of morphological components have been carried out and colours were assigned to localised saturated clusters, as shown Fig. 1(c). Fig. 1(b) and (c) show that crystallisation indeed exists in our observation. It is foreseen that this localised structure indeed correlates to the resulting saturation states. Further correlation analysis between the local crystallised region and saturation stage is being carried out. Liquid bridges which are observed as ring-like structures as shown in Fig. 1(d) can be of interest in future study. However, at the present stage, only saturated regions larger than ߨ mm 2 are taken into account.

Results and discussion
In order to quantify the magnitudes of drainage-injection curves, liquid saturation (Sr) is defined as where ܸ ௪ and ܸ are the total volume of water and voids, respectively. Fig. 2 shows the liquid saturation calculated from the pump in comparison to the liquid saturation obtained from image-based processing for the entire experimental region for the drainage-injection cycles with different capillary number, Ca. The upper branch of each hysteresis curve is the injection section and the lower one is the drainage section. Due to the removal of trapped liquid volume, usually existing as liquid bridges between grains that can only be observed during drainage processes, the accumulated saturated region of drainage process is expected to be smaller than the accumulated area in injection process. The impact of Ca on different drainage-injection curves is driven by the contribution of drainage rates. When Ca (or the drainage rate) is extremely small, i.e. Ca = 1.44 × 10 ିସ , due to the slow flow rate, more water is trapped in the form of liquid bridges, the blue curve shows a noticeable difference compared with other cases. As Ca is increased, the curves are overlaid with each other and shifted closer to the green dashed line (points with liquid saturation extracted from pump and image equal to each other), the liquid saturation of drainage-injection range stays constant at around 0.25-0.85.
The saturated clusters are statistically analysed in terms of geometrical properties. The largest saturated cluster in each image, corresponding to the saturated region behind the continuous liquid front, is excluded from the statistical analysis of cluster sizes. Here, we first look at the effective size of the saturated clusters, defined as ݀ = ඥ‫ܣ‬ , where ‫ܣ‬ denotes the area of the cluster (area of localised saturated region). By rearranging the experimental data into multiple sub lists by a given bin size (0.91 mm), the size distributions of saturated clusters for drainage under different Bo* are depicted in Fig. 3. The probability density tends to the peak value around 2-3 mm. The decrease of probability is then generated by the sizes of clusters in the log-linear scale. A lognormal distribution function was selected here to represent the experimental data.
The experimental probabilities of the area distributions are fitted using the lognormal distribution function as (see Fig.4) where ߤ and ߪ are the mean and the standard deviation of the logarithm of random variables, respectively. Fig. 5 shows the relationship between lognormal distribution parameters and generalised Bond number Bo*. Since for both subfigures in Fig. 5 the absolute value of correlation coefficient ‫ݎ‬ is larger than 0.65, a relatively strong linear relationship lays in between these two parameters and Bo*. It is found that the scale parameter ߤ is inversely correlated to Bo*. During the drainage stage, it has been observed experimentally that at a large positive generalised Bond number, the liquid front quickly grows to a plateau, with small amount of clusters left behind. As the generalised Bond number decreases, narrowed branched fingers appear. For large negative values of Bo*, some large sized saturated clusters are trapped in the non-wetting phase due to the instability induced by viscous forces.

Conclusions
In this study we experimentally evaluated the saturated cluster size distribution and its evolution during drainage processes of wetted granular materials. The distribution is found to follow a lognormal distribution. It is concluded that the scale parameter ߤ and shape parameter ߪ negatively and positively correlate with Bo*, respectively. Due to the fact that a strong linear dependence exists between lognormal distribution parameters and generalised Bond number, Bo*, it is possible to generate a function containing ߤ, ߪ and Bo* for describing the cluster distribution at various saturation levels. Knowing this, combined with the total surface energy calculated by the area of interfaces, it is possible to include this additional grain-scale information in constitutive modelling of unsaturated soils using both the degree of saturation and generalised Bond number, Bo*.