Erika Shiota

Toshifumi Mukunoki

Laurent Oxarango

Anne-Julie Tinet

Fabrice Golfier

Micro-and macro-scale water retention properties of granular soils: contribution of the X-Ray CT-based voxel percolation method

Keywords: hanging column experiment, marker-controlled watershed, representative element volume, voxel percolation method

Water retention in granular soils is a key mechanism for understanding transport processes in the vadose zone for various applications from agronomy to hydrological and environmental sciences. The macroscopic pattern of water entrapment is mainly driven by the pore-scale morphology and capillary and gravity forces. In the present study, the drainage water retention curve (WRC) was measured for three different granular materials using a miniaturised hanging column apparatus. The samples were scanned using X-ray micro-computed tomography during the experiment. A segmentation procedure was applied to identify air, water and solid phases in 3D at the pore-scale. A representative elementary volume analysis based on volume and surface properties validated the experimental setup size. A morphological approach, the voxel percolation method (VPM) was used to model the drainage experiment under the assumption of capillary-dominated quasi-static flow. At the macro-scale, the VPM showed a good capability to predict the WRC when compared with direct experimental measurements. An in-depth comparison with image data also revealed a satisfactory agreement concerning both the average volumetric distributions and the pore-scale local topology. Image voxelisation and the quasi-static assumption of VPM are likely to explain minor discrepancies observed at low suctions and for coarser materials.

Introduction

Water distribution and transfer in the vadose zone plays a major role in a wide range of applications including agronomy for plant growth response to water availability [START_REF] Couvreur | Dynamic aspects of soil water availability for isohydric plants: Ffocus on root hydraulic resistances[END_REF], hydrology with water budget and infiltration or runoff assessment [START_REF] Ciocca | Effects of the water retention curve on evaporation from arid soils[END_REF], environmental sciences for contaminant transport [START_REF] Šimůnek | Contaminant tTransport in the Uunsaturated Zzone: Ttheory and Mmodeling[END_REF] and, more recently, soil mechanics for soil bearing capacity [START_REF] Imhoff | Physical Qquality Iindicators and Mmechanical Bbehavior of Aagricultural Ssoils of Argentina[END_REF]). These vadose zone processes are found to be mainly driven by the water retention curve (WRC), which is a fundamental property of soils relating the average capillary pressure (or suction) with the volumetric water content (θ). The WRC could be seen as a macroscale effective property averaging the complexity of water and air distribution inside the pore space and the wettability of the medium (including surface tension and contact angle effects). Although criticised for filtering an important part of the local-scale information [START_REF] Mcclure | A geometric state function for two-fluid flow in porous media[END_REF], the WRC still constitutes the basis of the majority of continuum models describing water equilibrium and flow in the vadose zone.

Conventionally, the hanging column test (HCT) is conducted to evaluate the water retentivity of sandy soil. It consists of applying successive suction steps and measuring the change in sample mass due to the water drainage (or imbibition) at hydrostatic equilibrium. However, it is limited by the range of suction pressure that can be applied. Pressure cell or pressure plate methods can be an alternative for measurement at low water content but the accuracy of air pressure regulation is usually inappropriate for coarse materials. Moreover, the HCT requires very careful handling to guarantee appropriate reproducibility.

In order to overcome these issues, several methods have been proposed to indirectly estimate the WRC. Mercury intrusion porosimetry (MIP) provides an estimation of soil pore (or throat) size distribution (PSD). A straightforward use of Laplace's law is then commonly used to derive the WRC. Hence, pore structure of soil is generally modelled as a bundle of parallel pipes with sizes corresponding to the PSD. Such an approach thus fails to take into account the complex pore structure of the medium (e.g. tortuosity and constriction effect). To overcome these limitations, numerical derivations of the percolation theory have been proposed by [START_REF] Wildenschild | X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems[END_REF] and [START_REF] Bird | Simulation of soil water retention using random fractal networks[END_REF] in order to include, an estimation of the pore network complexity while deriving the WRC from MIP. However, the mathematical background of the percolation theory limits these approaches to 2D configurations.

In parallel, for more than a decade, the micro-focused X-ray computed tomography (µCT) scanner has proven to be a powerful tool to study soils at a micrometric resolution (see the recent reviews of [START_REF] Cnudde | High-resolution X-ray computed tomography in geosciences: aA review of the current technology and applications[END_REF][START_REF] Wildenschild | X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems[END_REF]. Therefore, this non-invasive imaging technique is able to extract the complete pore structure of sandy soils.

Moreover, µCT provides detailed information about the phase distribution and topology within pores that give insight not only on the fluid saturation values, but also on local properties such as phase connectivity, contact angle or interfacial curvatures (Al-Raoush and Willson 2005; [START_REF] Armstrong | Linking pore-scale interfacial curvature to columnscale capillary pressure[END_REF][START_REF] Andrew | Pore-scale contact angle measurements at reservoir conditions using X-ray microtomography[END_REF][START_REF] Armstrong | Beyond Darcy's law: tThe role of phase topology and ganglion dynamics for two-fluid flow[END_REF][START_REF] Manahiloh | Determining the sSoil Wwater Ccharacteristic Ccurve and Iinterfacial Ccontact Aangle from Mmicrostructural Aanalysis of X-rRay CT Iimages[END_REF]. [START_REF] Hamamoto | Pore network structure linked by X-ray CT to particle characteristics and transport parameters[END_REF] derived the mean pore diameter and pore coordination number by percolation thresholds and pore network connectivity from µCT observation on repacked columns of sand and glass beads representing different size fractions and particle shapes at different compaction levels. Recently, [START_REF] Manahiloh | Determining the sSoil Wwater Ccharacteristic Ccurve and Iinterfacial Ccontact Aangle from Mmicrostructural Aanalysis of X-rRay CT Iimages[END_REF] measured the WRC and values of interfacial contact angle on an unsaturated glass beads sample in a suction-controlled setup directly imaged with µCT.

These recent advances in µCT imaging, combined with theoretical or computational models, bridge the gap between pore-scale observations and macro-scale properties. Column-scale values of capillary pressure were back-calculated by [START_REF] Armstrong | Linking pore-scale interfacial curvature to columnscale capillary pressure[END_REF] from local measurement of interface curvatures during drainage and imbibition processes and comparison with direct measurement has shown good agreement. [START_REF] Berg | Connected pathway relative permeability from pore-scale imaging of imbibition[END_REF] carried out Navier-Stokes flow simulations on synchrotron-based µCT pore-scale images of an imbibition experiment in sandstone rock to determine relative permeability curves.

However, numerical computations of core scale parameters like retention curve or saturated and unsaturated hydraulic conductivity require to achieve a representative elementary volume (REV) of the soil sample, which is itself dependent on the targeted property [START_REF] Blunt | Pore-scale imaging and modeling[END_REF][START_REF] Raeini | Direct simulations of two-phase flow on micro-CT images of porous media and upscaling of pore-scale forces[END_REF]. [START_REF] Georgiadis | Pore-scale micro-computed-tomography imaging: nonwetting-phase cluster-size distribution during drainage and imbibition[END_REF] suggested that the two-phase REV for a sintered glass-bead porous medium is one to two orders of magnitude larger than the porosity-based REV and not reachable by µCT imaging at the required resolution. A similar conclusion was drawn by [START_REF] Yuan | Pore-scale simulations of drainage in granular materials: fFinite size effects and the representative elementary volume[END_REF] from pore network simulations on granular materials. In contrast, REV calculations by [START_REF] Costanza-Robinson | Representative elementary volume estimation for porosity, moisture saturation, and air-water interfacial areas in unsaturated porous media: dData quality implications[END_REF] on different unsaturated homogeneous porous media with varying degrees of heterogeneity in moisture distributions claim that moisture saturation-based and air-water interfacial area-based REVs are overly restrictive and requiring a porosity REV to be achieved can be sufficient. Nonetheless, two-phase flow simulations at the pore scale are still numerically costly and concerns remain regarding the capability to image soil samples at the REV size and directly compute the core-scale properties of these 3D digitalised porous media. In this framework, morphological approaches represent a promising alternative to exploit large volumes of 3D µCT images but at the cost of highly simplified physical description [START_REF] Hilpert | Pore-morphology-based simulation of drainage in totally wetting porous media[END_REF][START_REF] Vogel | Comparison of a lLattice--Boltzmann mModel, a Ffull-mMorphology mModel, and a Ppore Nnetwork Mmodel for determining capillary pressuresaturation relationships[END_REF][START_REF] Mukunoki | X-ray CT analysis of pore structure in sand[END_REF]. [START_REF] Berg | Connected pathway relative permeability from pore-scale imaging of imbibition[END_REF] showed that the morphological approach produced a good prediction of the relative permeability curve of sandstone rock for drainage steady-state experiment but it failed for imbibition conditions. Recently, [START_REF] Mukunoki | X-ray CT analysis of pore structure in sand[END_REF] introduced an original method called the voxel percolation method (VPM), which enables an a priori estimate of the retention curve based on the solid-pore binary image obtained from µCT. The WRC predicted by VPM compared fairly with experimental WRC, provided that the computation was carried out in a sufficiently large image sample.

These recent works suggest that in-depth assessment of morphological approach capabilities is required to gain confidence in outcomes of numerical models applied to soil drainage.

Notably, it remains unclear whether these approaches are limited to predicting WRC or whether lower-scale information can be properly captured in spite of their simplified mechanistic description. The present work aimed to address these issues to ultimately obtain a better understanding for the water retention property of granular materials.

To achieve this goal, three different materials were selected and their pore structure and water distribution were visualised by µCT during HCT retention curve characterisation. The VPM model [START_REF] Mukunoki | X-ray CT analysis of pore structure in sand[END_REF]) was considered as representative example of morphological approach. In the first part of this paper, the drainage WRC determination using a mini HCT apparatus is described. In the second part, a REV analysis is performed to establish the confidence associated with the properties of interest (porosity, water content and surface or interface specific areas) and the representativeness of samples with respect to saturation state. In the third part, the water-filled pore size distribution (WPSD) is analysed with respect to the suction state. Finally, the VPM is applied to estimate the WRC from the initial state binary image. This numerical model is then compared with experimental data both at the macro and pore scale.

Materials and methods

Experimental method

In order to evaluate the applicability of the method to different granular media, three different granular materials were selected. Toyoura sand is a clean silica sand commonly found in Japan.

It has a uniform fine grain distribution of median diameter 160 µm. Keisa5 sand is a silica sand from Japan with a uniform grain distribution of median diameter 500 µm, larger than Toyoura sand. Toyoura and Keisa5 sands also contain iron material. The third material is glass beads (CGB) manufactured by Potters-Ballotini, which were crushed in order to create an artificial material with very angular particles with a median diameter of 650 µm. Characteristics of the three materials are presented in Table 1.

A miniaturised HCT apparatus was designed to measure the WRC of the different materials.

The miniaturisation allows imaging of the whole sample with a maximum resolution of 6 µm with the X-ray micro-CT scanner. The sample of 15 mm-diameter was set up in a cylindrical cell, above a glass filter (20-30 µm pore diameter) (Fig. 1b). The sand was poured in the waterfilled cell to avoid air entrapment. The sample was compacted reaching the dry densities shown in Table 1. Suction (i.e. negative water head) was imposed by modifying the relative position of the top of the sample and the free water level in the syringe (Fig. 1a) in a range between 0 and 60 cm. It was measured based on a graduated ruler with an accuracy of 1 mm. It should be noted that selected suctions, aside from the initially saturated conditions, were significantly larger than the sample height. The distribution of suction steps for each sample was selected after the preliminary experiments in order to match the suction values that could be simulated with the VPM. The top of the sample was selected as reference to ensure a total saturation for a suction of 0 cm. The volume of drained water in the syringe was measured with a precision of 0.02 mL.

For each suction step, the equilibrium was guaranteed by a waiting time of at least 24 h. The experiment is performed under temperature-controlled conditions (20 ± 1°C).

For each suction step, once the equilibrium is reached, the sample was thermally insulated using polystyrene and hydraulically isolated by closing the outlet valve. The sample was then fitted in the µCT scanner (Toshiba TOSCANNER 32300 FPD) to perform the image acquisition around the centre of the sample. The resolution used was of 10 µm for crushed glass beads and Keisa5 sand, and 6 µm for Toyoura sand. Details on the µCT scan conditions are shown in Table 2.

Image treatment and analysis

All image treatment steps were performed on the open-source software Fiji [START_REF] Abramoff | Image pProcessing with ImageJ[END_REF][START_REF] Schneider | NIH Image to ImageJ: 25 years of image analysis[END_REF] or using custom-made programs. 3D grey-scale images were reconstructed from each set of 2D slices obtained using µCT. The grey value of each voxel depends on the material parameters such as density and atomic number, therefore allowing differentiation of the various materials present in the sample such as silica grains, water and air. However, when using the µCT scanner, artefacts may occur as well as image noise that influence the segmentation process [START_REF] Wildenschild | X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems[END_REF]. As a consequence, filtering of the images is necessary before the segmentation process. In the present study, fi 3D Gaussian blur with a standard deviation of 1 in order to remove noise and facilitate edge detection [START_REF] Tamura | Kkonpyuta gazoushori [Computer image processing[END_REF][START_REF] Wildenschild | X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems[END_REF].

Numerous segmentation approaches exist to distinguish different phase in an imaged porous medium [START_REF] Wildenschild | X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems[END_REF][START_REF] Berg | Generation of ground truth images to validate micro-CT image-processing pipelines[END_REF]. Recent works on multiphase porous medium have used region-growing segmentation [START_REF] Hashemi | A tomographic imagery segmentation methodology for three-phase geomaterials based on simultaneous region growing[END_REF][START_REF] Schlüter | Image processing of multiphase images obtained via X-ray microtomography: aA review[END_REF][START_REF] Mukunoki | X-ray CT analysis of pore structure in sand[END_REF]). In the present work, a custom-made marker-controlled watershed software [START_REF] Soille | Morphological image analysis: Pprinciples and Aapplications[END_REF][START_REF] Mukunoki | X-ray CT analysis of pore structure in sand[END_REF] was used. This segmentation technique has already shown a good capability to separate three phases in unsaturated sand samples while recovering volume fractions matching gravimetric measurements [START_REF] Mukunoki | X-ray CT analysis of pore structure in sand[END_REF].

Marker-controlled watershed is based on seeds for each phase and iteratively simulates their growth so as to fill the image (Fig. 2). The resulting discrete image was analysed considering three phases: water, air and solid (Fig. 2). In Toyoura and Keisa5 sands, the solid phase consisted of the sum of silica and iron grains. The sum of water and air phases is referred to as the pore phase. Segmented images were checked to validate the method based on visual inspection of 2D images with superposition of pre-and post-treated data once every 50 images, especially regarding separation between the pore and solid phases.

For the subsequent studies, the final image was of the cylindrical cell borders, which simplifies the numerical treatment as well as providing improved computing performances. Morphological characteristics were computed in order to validate the image processing through comparison of image data and experimental results, but also to quantitatively study the evolution of unsaturated properties depending on the saturation level and material properties.

Global characteristics under consideration include volume averaged properties (porosity (ε), θ , specific surface and specific air-water interface) and statistical description of local geometrical properties (WPSD and air-filled particle size distribution (APSD)) which allows a more detailed description of the sample structure.

The ε and θ were straightforwardly calculated from the segmented images as the ratio between the number of corresponding voxels (respectively pore and water voxels) to the total number of voxels in the sample. Interfacial areas are evaluated through the discrete cutting lines method proposed as a 3D extension of Cauchy-Crofton theory [START_REF] Santalo | Integral geometry and geometric probability[END_REF][START_REF] Legland | Computation of Minkowski measures on 2D and 3D binary images[END_REF][START_REF] Lehmann | Efficient N-Ddimensional surface estimation using Crofton formula and run-length encoding[END_REF]. This method is less sensitive to the staircase effect than a direct counting of interface voxels and has much lower computational cost than methods based on surface meshes. The specific surface (Ss in m-1) was computed as the ratio between the area of pore-solid interfaces and the total volume of the sample. In a similar manner, the specific water-air surface (Swa in m-1) was computed using the area of water-air interfaces.

Pore size distributions (PSD) were calculated using the 3D granulometric method (GM) developed by [START_REF] Mukunoki | X-ray CT analysis of pore structure in sand[END_REF]. This method considers pores as spheres of various radii.

For each radius, the occupation rate of spheres having three points of contact with the solid phase is computed. PSD is derived from this occupation rate and the diameters of the spheres. In a classical manner, GM is used to estimate the total PSD (TPSD) of the porous medium. In this study, it was also applied to the water phase only to estimate the WPSD.

A major concern with regards to studies at the pore structure scale is the representativeness of the selected medium in relation to a larger scale, relevant to the applications. To consider this issue and allow upscaling, the evaluation of the REV is necessary. In order to estimate a representative property of the medium, the study should involve a volume large enough to capture the average trends at the upper scale. All the macroscale models used to describe soils (and more generally porous media) can be used confidently only if a REV can be defined for the studied process. Various methods have been proposed to estimate the size of the REV using a reconstructed 3D image of the pore structure. Of these, the stochastic REV [START_REF] Zhang | Pore scale study of flow in porous media: scale dependency, REV and statistical REV[END_REF] was selected in the present work. This method offers the advantage to provide a distribution of the studied property for each region of interest (ROI) volume. The average value of the distribution generally converges towards the representative value, whereas the discrepancy of the distribution, computed as the relative standard deviation (RSD), provides a metric of the representativeness of the studied ROI size. In the present study, distributions were computed as a set of 20 ROI distributed stochastically in the total volume. This ensured that the discrepancy between the total volume property and the distributions means was lower than 1%. The ROI As emphasised in the Introduction, REV is not an intrinsic property of the medium. Its size generally varies depending on the studied property. Therefore, REV should be determined for each studied property of the medium. In the current work, calculation of the REV was carried out for each global characteristic (porosity, water content, specific surface and air-water interface).

Numerical modelling

The numerical derivation of the drainage WRC is achieved using the VPM model [START_REF] Mukunoki | X-ray CT analysis of pore structure in sand[END_REF]. It is based on two steps. The first, using the same algorithm as the 3D GM, consists of labelling each pore, considered as a sphere element, with the sphere geometrical characteristics. Once this cluster labelling process is carried out, a percolation step can be performed for different values of the capillary pressure. The Young-Laplace law is used to associate a pore diameter to this capillary pressure with the water surface tension at 20°C (0.07288 N m-1) and a contact angle for the air-water-quartz system of 44° ± 5° (based on data from Ethington 1990;[START_REF] Mukunoki | X-ray CT analysis of pore structure in sand[END_REF]. During each percolation step, fluid is allowed to penetrate, from the entry face into pores labelled with spheres of diameters higher or equal to the chosen value as long as the injected fluid phase remains connected. Fluid-filled voxels may then be marked so as to allow calculation of the injected fluid content. If this calculation is performed decreasing the size of the sphere diameters, the drainage WRC is derived.

Conversely, if the calculation is carried out with increasing spheres, the imbibition WRC is derived. It should be noted that the capillary pressure points that can be computed with VPM are strongly dependent on the image resolution since diameter values can only be multiples of the voxel size.

In the present work, VPM simulations were performed based on the initial state (saturated conditions) image for each sand. The confidence interval associated to each VP simulation was based on the range of contact angle variation (between 39° and 49°), with a larger contact angle leading to a smaller air entry pressure of the WRC. As for µCT image analysis, PSD calculations were performed on VPM simulations for comparison purposes. Note that only the total and APSD were directly computed here. The WPSD was obtained as the difference between TPSD and the APSD.

Results and discussion

Macro-scale results: the WRC

The experimental WRC (based on water outflow measurement) for the three materials is presented in Fig. 3 along with the image-based WRC (based on voxel counting), and the WRC predicted by the VPM. All WRC are typical of sandy materials with low air entry pressures, which decrease with increasing mean grain sizes. A strong decrease of θ was then observed at intermediate suctions leading to small residual water content (< 0.05). The variability of the VPM data with contact angle decreased with suction and, thus, with the increase of water content.

Experimental and image-based WRC show very similar global behaviour, those values are mostly included within the range provided using VPM. Some differences should, however, be noted. Since the capillary pressure is not recalculated from the image, water content between imaged WRC and experimental WRC may be compared at each drainage step. Whilst the comparison for Toyoura sand is rather good, it is less satisfactory for Keisa5 and CGB. Indeed, variation in water content of around 0.05 was observed for Keisa5 Step 2 to 4 and CGB Step 3 and 4. We believe that the sub-sampling procedure could explain this bias due to preferential pathways for air intrusion along the walls of the experimental sample, especially for the coarse materials.

The VPM data tend to deviate from experimental data at low water content. Indeed, under such conditions, air invades the smaller pores, making the VPM much more sensitive to the image resolution and potentially leading to a larger calculation error.

Evaluation of the REV on averaged characteristics

In this section, we investigate the changes in the REV size achieved when the RSD is below 5%. Since there is no consensus for the best value of RSD to describe the REV, the value was arbitrarily chosen. Besides, the ROI size was limited to 350 voxels (i.e. half of the total subsample size in each direction). The evolution of RSD with ROI size for each sample at each step considering ε (initial step of volume REV), θ (step 1 to 4 of volume REV), specific surface Ss (initial step of surface REV) and Swa (Steps 1 to 4 of surface REV) are shown in Fig. 4.

Total values of each studied characteristic as well as computed REV and RSD for the maximal ROI size are summarised in Table 3.

The computed Ss ranging from 6661 m-1 for CGB up to 16817 m-1 for Toyoura sand appears to scale properly with the median grain diameters and the median pore size for each material. Contrasted behaviours were observed for each material concerning Swa. Except the null value at Step 1, reflecting the full saturation, Toyoura sand exhibited a constant Swa of ~2300 m-1 for Steps 2 to 4, while θ is continuously decreasing. This reflects a transition from a small number of large water clusters to a large number of small clusters. For Keisa5 sand, the same effect appeared to promote a maximum Swa value at the intermediate water content of

Step 3 (14.4%). The very low water content at Step 4 (4.6%), leading to the reduction of the number of water clusters, probably explains this noticeable decrease of Swa. Finally, CGB exhibited very low Swa from Step 1 to 3, which is in agreement with its high water contents. At

Step 4, Swa reached a value of 731 m-1, which is consistent with the results obtained for Keisa5 sand due to the similar mean grain size for both materials.

RSD evolutions with ROI size for volume (Fig. 4a,c,e) and surface (Fig. 4b,d,f) properties yield very similar trends with the notable exception of samples close to water saturation. In these cases, the air-water specific interface is limited due to the small amount of air. Therefore, the likeliness to have some air phase in a given ROI is small, leading to a very large RSD for Swa (over 100%).

In saturated conditions (initial step) or at high water saturation, both REV values are reasonably low for all materials (i.e. below 4 times the grain size, Table 3). This value, however, increases with the progress of the drainage process to the point that neither volume REV nor surface REV with RSD lower than 5% is achieved at Step 4, with the exception of CGB surface REV. However, in the latter case, the water content is much larger than for Keisa5 at Step 4 and Toyoura at Step 4 (13.5% against 4.6% and 7.1% respectively). The REV size globally evolves in a similar way as the correlation length of the object of interest (i.e. the water clusters). As the water clusters become increasingly fragmented when θ decreases, the inter-cluster characteristic length increases and hence the REV increase. Note, however, that at high water content, the estimate of this correlation length may be biased because of the cut off effect of the experimental sample size [START_REF] Georgiadis | Pore-scale micro-computed-tomography imaging: nonwetting-phase cluster-size distribution during drainage and imbibition[END_REF].

Considering the value of RSD for ROIs of 350 voxels generally below 15% (Table 3) and the decreasing trend observed in Fig. 4, it is reasonable to consider that the ROI of size 600 x 600 x 600 and 700 x 700 x 700 voxels used for characterisation purposes present satisfactory representativeness for the three materials. These ROI sizes are 3-7 times larger than the porosity REV size, which was considered large enough to study multiphase properties by [START_REF] Costanza-Robinson | Representative elementary volume estimation for porosity, moisture saturation, and air-water interfacial areas in unsaturated porous media: dData quality implications[END_REF]. In the present study, RSD values measured for the porosity REV were generally larger than 30% for all the computed properties and thus appear inadequate for a confident interpretation of experimental and numerical results. However, it must be kept in mind that defining a REV for multiphase property presents an intrinsic difficulty due to the increasing size of the connected phase with increasing sample size, as suggested by [START_REF] Georgiadis | Pore-scale micro-computed-tomography imaging: nonwetting-phase cluster-size distribution during drainage and imbibition[END_REF] and [START_REF] Yuan | Pore-scale simulations of drainage in granular materials: fFinite size effects and the representative elementary volume[END_REF].

Pore (TPSD) and water-filled (WPSD) pore size distributions

The PSD values were calculated at the different steps for each material from the segmented images as well as from VPM simulations and the comparison is shown in Fig. 5.

Under initial conditions, the WPSD corresponds to the TPSD. Moreover, since all VPM simulations were performed on the initial images, TPSD from image data and VPM data were strictly identical. All materials show a unimodal distribution with increasing dispersion around the mode from Toyoura to Keisa5 and to CGB. More specifically, Toyoura PSD shows a mode of 66 µm pore diameter (110 µm for Keisa5 and 170 µm for CGB) with pores ranging from 5 to 150 µm (5-250 µm for Keisa5 and 30-350 µm for CGB).

For Toyoura sand (Fig. 5a,b), the WPSD exhibited very similar behaviours between imagebased and VPM-based WPSD. The successive suction steps led to a shifting of the right part of the WPSD since the largest pores were the first filled with air. The absence of a significant tailing effect confirms the good uniformity of the Toyoura sand and is consistent with its narrow distribution around the mean pore size. The VPM-based WPSD, however, exhibited a slightly stiffer shape, close to the Young-Laplace diameter used to define the air percolation at a given suction step.

For Keisa5 sand (Fig. 5c,d) and CGB (Fig. 5e,f) the tailing effect became increasingly marked on image-based WPSD with increasing air intrusion. At Step 3 of Keisa5 sand, some larger pores remained filled with water, even though smaller size pores began to be desaturated.

A similar observation was made on Step 4 of CGB for the imaged-based WPSD (Fig. 5 e,f).

The pore heterogeneity features of both materials were probably driving this effect; some large pores were isolated from the air flow by smaller size pores which played the role of throats or constrictions. Therefore, those large pores may only get filled by air when the capillary pressure is sufficient to invade these throats. The VPM only partially captured this behaviour, as can be water.

All in all, keeping in mind that the water contents are not strictly equal between the steps of VPM simulations and direct imaging, the VPM appears to satisfactorily describe the distribution of water within the pore system. However, being based on a straightforward use of the Young-Laplace diameter to impose a chosen suction, it tended to impose an overestimated cutoff of the WPSD close to this diameter. The resulting curves were thus steeper than those directly obtained from the µCT image. For material presenting a spread PSD, VPM simulations could partially describe the complexity of pore connectivity and the associated entrapment of water in large pores. However, dynamical processes, such as by-passing entrapment, preferential flow effects or unstable flow patterns, remain beyond the capability of the model and could explain some of the discrepancy between simulation and imaging [START_REF] Berg | Connected pathway relative permeability from pore-scale imaging of imbibition[END_REF].

Local phase distribution comparison

In order to go a step further in the analysis of the capability and limitation of the VPM, a straightforward comparison of the 3D phase fields (air, water and solid) obtained from direct imaging (I) and VPM simulations was performed. Considering similar water content conditions, the two datasets were compared voxel by voxel for Toyoura Steps 3 and 4, Keisa5 Steps 2-4 and CGB Step 4. 3D matching images were computed by defining nine different classes of voxels, each one consisting of the matching of the phases (air, water and solid) between the two sets of images (I and VPM). Fig. 6 details the proportion of each matching set expressed as the volume fractions of matching and mismatching voxels with respect to the total volume. The solid phase was mostly well matched, but some variation still remained between the segmented image and the VPM simulation ranging from 4.5% for Keisa5 Step 2 and 3 up to 8.5% for Keisa5 Step 4. Indeed, the VPM simulation was performed using the initial step binary image. Since the experiment was performed outside the µCT, the sample needed to be repositioned in the µCT for each image acquisition. Even though this repositioning was performed very carefully, considering the high resolution of the acquisition, even a very small modification of the initial location (e.g. tilting, rotation and modification of the structure) leads to some errors in the solid phase reconstruction. Moreover, it is well known that even if identical µCT acquisition is used, changes in the CTvalues could occur, resulting in a variability of the segmentation procedure. Such uncertainties promote the observed mismatches in the localisation of the solid phase.

With regards to the fluid phases, the dependency of matching success rate with the phase saturation is shown in Fig. 7 to ease the interpretation of Fig. 6. Note that all the comparison results are summed up in Fig. 7, independently of the material or of the saturation step. It may be observed that as the saturation of a given phase increased, the success rate also increased until reaching a plateau of 80% success rate from a saturation of around 50%. The remaining 20% could be caused by the grim resolution limit [START_REF] Saxena | Effect of image segmentaition & voxel size on micro-CT computed effective transport & elastic properties[END_REF] or the potential issue for segmentation [START_REF] Leu | Fast X-ray mMicro--Ttomography of mMultiphase Fflow in Berea sSandstone: aA sensitivity sStudy on Iimage Pprocessing[END_REF]. This maximum value of success rate was, for a large part, due to pore-solid mismatches, as discussed previously, but air-water mismatched voxels still remained. Such voxels represent between 4% and 9%, the latter value being for CGB Step 4, where the discrepancy in water content between segmented image and VP simulation was the largest.

The local distribution of matching and mismatching voxels in the 2D middle horizontal slice of the 600 x 600 x 600 subsample is illustrated in Fig. 8. A focus on 3D air voxels distribution is presented in 3D in Fig. 9.

The images of Toyoura sand at Steps 3 and 4 and Keisa5 sand at Step 4 (Fig. 8,9a,b,e) correspond to low water contents (0.11, 0.07 and 0.05 respectively). The agreement for the air phase, which is dominant, appears very satisfactory, as suggested previously by Fig. 7 but a reasonable agreement could also be noticed concerning the water phase. It was observed that the errors were rather homogeneously distributed through the volumes. Some of these errors may be attributed to the same causes than the pore-solid errors. VPM simulations appear to promote water films along solids and water bridges between two or more solids. Such geometrical anomalies may be attributed to imaging-induced roughness or voxelisation. Indeed, VPM considers a spherical shape for the air-filled pores, therefore voxels at the surface of the solids not matching the spheres or at junction between two spheres are not filled with air. It should also be noted that most of the pores anomalously filled with air by the VPM present a large size, whereas the small size pores filled with air in the segmented image remain occupied by water in the VPM simulation. This local observation supports the interpretation of the WPSD that strict application of the Young-Laplace theory prevents the invasion of pores smaller than the capillary diameter in VPM simulations. In practice, several phenomena occurring during the drainage experiment are likely to break this assumption. For instance, fluid flow is by essence associated with local pressure non-uniformity in both the water and gas phases during a drainage step. It is thus possible to locally reach capillary pressure higher than the imposed suction, thus allowing small pore invasion. Aside from these fundamental considerations, Fig. 9a, b puts forward a good uniformity of the satisfactory air-phase matching all along the sample height.

For mid-wet conditions (Keisa5 Step 3 and CGB Step 4), it was first observed that the water and air phases were rather heterogeneously distributed within the medium. This supports the hypothesis of preferential flow occurring during the experiment. Moreover, the air phase overestimation by VPM simulation was mostly attributed to mismatches with water phase locations, possibly due to the snap-off or bypassing mechanisms.

Finally, the 3D image for Keisa5 Step 2, in almost saturated conditions, showed the largest discrepancy with the μCT image (Fig. 8c and more drastically in 3D in Fig. 9c). The segmented image data exhibited non-connected air phase within the subsample. This indicates that the subsampling misses a large part of the air-phase intrusion. This is likely to be the main explanation for significantly lower experimental water content than imaged water content. The discrepancy between the segmented image and VPM air distribution may be imputed to boundary effect. Indeed, in VPM simulations, air is injected homogeneously over the whole entry surface (Fig. 9c), but in the experiment, the air distribution around the sub-sample image was dependent on the flow surrounding the selected ROI. It should also be noted that a rather low suction was applied to this sample (17.8 cm). The dominant effect of capillary forces considered in the VPM model may be questioned: gravity and viscosity effects could turn dominant and enhance air phase fingering.

Conclusion

In the present study, WRC were derived for three different granular materials. The samples were monitored using µCT at various steps of the drainage process during HCTs. A markercontrolled watershed method was applied to segment the different phases of each sample. Based on this segmentation, a morphological approach (VPM) was applied to simulate air intrusion on a subsample. Then experimental, imaged and simulated results were compared from the global (WRC and PSD) to the local scale (2D and 3D images of the two-fluid phases), including a REV analysis. The following results were obtained:

Fig. 1 .

 1 Fig. 1. Experimental set-up. (a) hanging column test experiment, (b) schematic of the cell and (c) set-up

Fig. 2 .

 2 Fig. 2. Marker-controlled watershed method. Top: computed tomography (CT)-image pre-processing of

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Water retention curve (WRC) for each sand (Toyoura, Keisa5 and glass beads (CGB)) using

Fig. 6 .

 6 Fig. 6. Matching phase volume against total volume between segmented image (iron phase is included 595

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Matching success rate between segmented image data and voxel percolation method simulation

Table 1 . Material properties in Toyoura sand, Keisa5 and CGB

 1

		Toyoura	Keisa5	CGB
	Median grain diameter (µm)	160	500	650
	Uniformity coefficient	1.5	1.5	1.5
	Dry density (kg L -1 )	1.51	1.72	1.53
	Initial porosity	0.43	0.34	0.39
	Grain density (kg L -1 )	2.64	2.63	2.51
	Sample Volume (mm 3 )	27.3	23.5	26.4

Acknowledgements

This study is part of the PROTINUS project funded by the European Union's Horizon 2020 research and innovation program under grant agreement No 645717. It was supported by the International Research Organisation for Advanced Science Technology (IROAST) and Priority Organisation for Innovation and Publisher: CSIRO; Journal: SR:Soil Research DOI: 10.1071/SR18179

The segmentation method showed good adequacy with the core-scale experimental measurements; however, there was limitation caused by the sub-sampling and the nonconsideration of preferential pathways for air intrusion along the walls of the test cell, especially for the coarser materials. The REV analysis supports the choice of subsample size (600  600  600 and 700  700  700 voxels) for the characterisation study of the three materials. Such ROI sizes are 3-7 times larger than the porosity REV size and are very likely to be of the same order of magnitude as the water content and surface REV, even if we have to be cautious with the definition of the REV in multiphase conditions due to the relative water cluster size compared with the core dimensions.

Through the above validation of image processing, this study proposes an original comparison of local phase field at the pore scale between direct imaging and VPM. Even though it is based on a simplified quasi-static approach. While considering the phase presenting the higher saturation (i.e. water at the beginning of the drainage or air at the end), the VPM provided a satisfactory agreement for the Toyoura sand sample. The agreement with direct imaging results tended to decrease for materials with larger PSD (Keisa5 sand and CGB). The assumption of capillary driven quasi-static drainage also promoted an overestimation of the cut-off effect close to the Young-Laplace pore diameter. Hence, the VPM model demonstrates a good capacity to reproduce air intrusion for the studied materials at the macro-scale.

Conflicts of interest

The authors declare no conflicts of interest.