Determination of the relationship among capillary pressure , saturation and interfacial area : a pore unit assembly approach

Three state variables namely, degree of saturation (S), capillary pressure (Pc) and specific air-water interfacial area (a) are indispensable for modelling coupled processes relevant to unsaturated soils mechanics, agriculture, and contaminant hydrology. They play a key role in simulating various phenomena and the determination of various parameters and physical characteristics such as the unsaturated soil shear strength, field capacity, wilting point, air and water diffusivity and the rate of dissolution of contaminants. The determination of soil water retention curve (S-Pc) as well as the specific interfacial area (a) using available experimental techniques is a challenging and time consuming task. Therefore, a numerical technique that employs basic soil properties to obtain these variables is of much value and high practical and theoretical importance. In the current study, the porous network extracted from a discrete element model (the so-called pore unit-assembly) has been used to directly model the drying and wetting processes inside a granular soil packing and to obtain the values of Pc, S and a. The results of the simulations are in good agreement with the experimental data, which points to the efficacy and adequacy of the introduced algorithms and involved assumptions for this purpose.


Introduction
The soil water retention curve is a fundamental relationship in unsaturated soil mechanics, based on which different hydraulic and mechanical parameters are defined.Knowledge of this curve is essential for determining the contribution of matric suction to the shear strength in unsaturated soils and obtaining different hydraulic parameters such as water diffusivity.It is a key part of models used for simulating water flow and solute transport as well as soil mechanical behaviour [1].
The soil water retention curve, however, is not unique and for different hydraulic paths (drying, wetting, and scanning), different equations need to be introduced to capture variation of capillary pressure versus degree of saturation.This phenomenon has been referred to as hydraulic hysteresis [2][3][4][5][6].There is also hysteresis in mechanical properties such as effective stress parameter [7].
The hydraulic (and associated mechanical) hysteresis can be attributed to the fact that at different hydraulic paths and the same degree of saturation, the water is distributed differently inside the medium [8].In other words, at the same saturation, different macroscopically measured capillary pressure values are found due to different distributions of water phase (manifested in different air-water interface configurations).Therefore, in addition to saturation, a macroscale variable is required to account for many possible distributions of water at a given saturation.In other words, a new state variable is needed for hydro-mechanical models.This is proposed to be the specific interfacial area, defined as the amount of air-water interfacial area per total volume of an REV.The inclusion of this variable in coupled hydro-mechanical models can allow for modeling hysteresis.That is, the hysteretic capillary pressure-saturation curves should be replaced by a unique capillary pressure-saturationinterfacial area surface [7][8][9][10].
The measurement of capillary pressure and saturation is already a time consuming task.For the investigation of the relationship among capillary pressure, saturation and interfacial area, one needs also to deal with the measurement of specific interfacial area.It adds further challenge and difficulty and is often done through sophisticated experimental methods such as tracer technique (e.g.[11]) and x-ray micro-tomography (e.g.[12]).Another approach for the study of role of interfacial area has been the use of pore-scale computational models such as lattice Boltzmann [13,14], SPH [15], and pore-network modelling [16].These models, however, often suffer from serious difficulties to be used for practical purposes in geotechnical analysis.They need either the soil porous network or are highly computationally demanding, and moreover, they cannot easily and directly make use of the soil basic properties such as porosity and grain size distribution.
The purpose of this study is to present a method to use a grain scale model (DEM) rather than a pore-based model.The merit of such an approach is that the grain size distribution and porosity can be directly used when constructing the model of the granular packing and then the porous network would be obtained from the packing.Next, wetting and drying are simulated therein and the required state variables (air-water interfacial area, capillary pressure and saturation) are obtained.
In what follows, the required steps, namely, packing construction, extraction of porous network, and invasion criteria for wetting and drying are described and simulations versus experimental data are presented.Finally, concluding remarks and future perspectives are presented.

Model description 2.1 Model assumptions and considerations
In order to start with a medium where less microstructural and mineralogical complexities and heterogeneities are present to affect a aw -S w -P c relationship, this study is focused on modelling a packing consisted of glass beads.Therefore, at this stage, our focus is mainly on the coarse-grained material.We consider rigid grains and disregard possible mineralogical effects (such as water absorption by the grains and reaction with the grains).Other assumptions are that grains are modelled as spheres and the contact angle is assumed to be equal to 0° for drying and wetting.The contact angle hysteresis will be added in future investigations.

Construction of packing
For constructing the packing, we use the open source discrete element model YADE-DEM [17].Discrete element technique allows for adjusting particle packing while the grains are allowed to move and interact.We can thus start with grains at arbitrary locations and let them move and reach a stable configuration under their own weight as well as external pressure.We take the model domain to be a box, with dimensions 0.02m x 0.02m x 0.02m, within which we randomly position a cloud of 4000 particles.Grain sizes are chosen from a particle size distribution for the glass bead packing used in drying-wetting experiments of [18].Values of mechanical parameters of the particles used in our simulations are reported in Table 1.The values of mechanical parameters of the domain boundaries were the same as those for the particles.In this study, Hertz-Mindlin contact law has been used.
The initial porosity of the cloud of particles is larger than 0.99.To compact the cloud of particles into a loose packing, we apply a confining pressure of one kPa on all boundaries of the box.The model is then run until the cloud of particles is compacted and a force equilibrium is reached, which corresponds to an unbalanced force smaller than 10 -6 N. The packing in equilibrium has a porosity value of 0.34.The next step is to simulate an oedometer test.For this purpose, we fix all boundaries except for the top boundary.A confining stress of 1MPa is applied to the top boundary to simulate compaction of the packing.Further details on the utilized discrete element approach can be found in [17].

Table 1.
Input micromechanics parameters for glass beads (Chung and Ooi, 2008) [19] * Different value selected to match the porosity value of 0.34 reported by [18]

Extraction of pore network of the packing
The pore structure of the packing needs to be extracted in order to simulate the wetting and drying processes.The pore network of a granular packing consists of pore bodies and throats.Pore bodies constitute the greatest portion of total void volume.They are connected by throats which are considered to be without any volume (either as 2d surface elements or 3d but with negligible volume).There are various pore network extraction techniques originally introduced for extraction of the pore network from three-dimensional images of porous media [20].Here we use a regular tessellation technique (see [21] and [22]).The triangulation results in the division of pore space into several tetrahedra, each containing one pore unit.Vertices of each tetrahedron are located at the centre of the four particles embracing a pore unit (Figure 1).The scalene triangle formed between two neighbouring tetrahedron is considered as the pore throat connecting the two.

Invasion algorithms
Following Yuan et al. [23], we have used the Mayer-Stowe-Princen (MSP) [24][25][26] to determine throat entry pressure surpassing which the air invades a pore throat and connected pore units during drainage.For imbibition, we assume that the critical curvature at which water fills a pore unit is that of the inscribed sphere of the pore unit.
The capillary pressure corresponding to the radius R i of the inscribed sphere in pore unit is given by: where γ is the interfacial tension between water and air (72 dynes/cm 2 ).However, we need to account for cooperative filling during imbibition; the more the number of air-saturated neighbouring pores, the smaller the capillary pressure at which imbibition occurs.To account for this, we follow the method of Jerauld and Salter (1990) [27] to modify equation ( 1): where N air is the number of air-saturated neighbouring pore units.
The radius of the inscribed sphere is found by using the Cayley-Menger determinant, which is a method for finding the radius of an inscribed sphere for four touching or non-touching particles.For more information, the reader is referred to MacKay (1973) and Michelucci and Foufou (2004) [28,29].

Interfacial area determination
To determine the air-water interfacial area, we calculate the area of main menisci, also-called "terminal menisci", of interfaces inside a throat and pore unit.We neglect the water films on the solid surfaces and the liquid bridges that may be present between two touching particles.The position of an air-water interface (and thus, the amount of interfacial area) depends on the capillary pressure.Let's assume that an interface moves from pore unit i towards pore throat ij during drainage, and the reverse occurs during imbibition.Thus, capillary pressure varies between a minimum value of P ci =2γ/R i and a maximum value of P cij =γ/R ij.Then, the area of the meniscus varies between a maximum and a minimum value at these locations, respectively.The minimum value of the the air-water interfacial area is simply equal to the area of the scalene triangle that forms the pore throat ij.The maximum value of the the air-water interfacial area is simply equal to the area of the meniscus formed on the inscribed sphere of pore body i.It has an area ‫ܣ(‬ @ ௪ ) defined by: where ܸ ∩ is the volume of a portion of pore unit i that is assigned to pore throat ij and V i is the volume of pore unit i.In between these two extremes, we assume that the air-water interfacial area changes proportional to (γ/P c ) 2 .During imbibition, it may happen that an interface remains hanging on the inscribed sphere of a pore unit, as the inscribed sphere has a much larger radius than the pore throat.When the capillary pressure is smaller than the entry pressure of pore throat ij but larger than the capillary pressure corresponding to the inscribed sphere, then we assume that the air-water interface remains hanging on the narrowest opening of pore throat ij.When the capillary pressure is lower than that of the inscribed sphere, spontaneous imbibition occurs which fills the pore unit.This will create new interfaces in other pore throats.Finally, the interfaces that are located next to a residual phase are assumed to be located exactly at the narrowest opening of a pore throat

Results and discussion
Figure 2 shows the pore unit and throat size distributions of the pore unit assembly extracted from the modelled packing.Figure 3 illustrates the capillary pressuresaturation curve of sample with porosity of 0.34 as compared to the experimental data.The agreement between the simulated and experimental data is not fully satisfactory.Although in most parts of the drying and imbibition curves a good agreement is seen, the experimental residual saturation is not fully captured.This can be attributed to phenomena which have not been accounted for in our model such as film flow or the presence of liquid bridges, and their coalescence.
The variation of air-water specific interfacial area versus saturation has been depicted in Figure 4, where experimental data are also plotted.Here also, we find discrepancy in estimating values of air-water interface at low degrees of saturations.
Furthermore, the variation of air-solid specific interfacial area, a as , versus saturation are presented in Figure 5; as it can be seen, less hysteresis is present in this relationship as compared to the air-water specific interfacial area versus degree of saturation.a as versus saturation shows almost a linear relationship with degree of saturation.This is in agreement with the results of Culligan et al. (2004) [18].

Concluding remarks
In this study, a method for the determination of interfacial areas in an assembly of pore units and throats was presented.The porous network was extracted from granular packing by regular tessellation.Our results were of reasonable accuracy as compared to the experimental data.However, further comparison with experimental data of other packings and granular materials (e.g.sand) and additional validation is necessary to assess the scope of application and the accuracy of the proposed method.Inclusion of further complexities to the model such as including contact angle hysteresis, as well as required local rules (relationships) to account for film flow, trapping and snap off would enhance the capabilities of the model for a better estimation of the values of airwater interfacial area and its variation with capillary pressure and saturation.

Figure 1 .
Figure 1.(a) One pore unit embraced by four particles; (b) Inscribed circle in a pore throat, with radius R ij (c, d) Inscribed sphere in a pore body with radius R i .

Figure 2 .
Figure 2. Pore throat radii (Rij) distribution and the distribution of inscribed spheres in pore units (Ri).

Figure 5 .
Figure 5. Simulated air-solid specific interfacial area versus saturation