Drying-induced volumetric behaviour of clays interpreted via binary pore-scale modelling

The volumetric response upon drying is relevant in the geotechnical and ceramic industry due to the key role played in crack formation. The volumetric response of clays upon drying is relatively complex; the void ratio decreases considerably in the saturated state, levels off almost abruptly at the onset of desaturation, and increases at very low water contents. This paper presents an original deformable pore-network model (DPNM) to elucidate the micromechanisms controlling the complex volumetric response of clays upon drying. The main features of the DPNM are the ‘binary ’ approach (pores in the clay are assumed to be either fully saturated with water or dry) and the use of independent tests to calibrate the parameters of the model (compression behaviour of saturated reconstituted clay and dry powder). The DPNM shows that all pores shrink as pore-water tension increases in the saturated state, although the larger pores contribute more to the overall reduction in void ratio. At the onset of desaturation, the larger pores dry out and rebound, whereas the smaller pores remain saturated and shrink. At very low water content/high pore-water tension, the desaturation of the smaller pores causes them to rebound, thus generating a macroscopic increase in the void ratio.


Introduction
The hydromechanical behaviour of clay upon drying is a complex process, as it spans a wide range of degrees of saturation and water contents (from plastic to limited shrinkage).Experimental observations of clays dried from the slurry state to very high pore-water tension (suction) show a complex volumetric behaviour.The void ratio decreases considerably in the saturated state, levels off almost abruptly beyond the desaturation point (at degrees of saturation above ~90-95 %, TARANTINO et al., 2010), and increases (volumetric rebound) at very low water contents (Tarı ̀ and Ferreira, 1998, Tarı ̀ et al., 1999, Cunningham et al., 2003, Vesga, 2009).To the author's best knowledge, there are no conceptual/numerical models encompassing all these aspects of the volumetric response of clays upon drying.
Volume deformation and deformation-controlled boundary conditions are directly associated with the formation of cracks (Cordero et al., 2017), which can compromise the performance of geotechnical structures due to their effects on soil mechanical and hydraulic behaviour.The presence of tension cracks influences the stability of natural slopes and flood embankments due to water seeping into the desiccated crest (MARSLAND, 1971;).Shrinkage behaviour is also important in the ceramic industry, as cracks develop in wares during the drying of slip-cast forms (; Tarantino et al., 2010).
This paper presents an original Deformable Pore-Network Model (DPNM) focused on unconstrained drying to elucidate the micromechanisms controlling the complex volumetric response of clays upon drying.PNMs have been shown to be a reliable tool for investigating multiphase flow in porous media (Fatt, 1956, Joekar Niasar et al., 2009, Rajaram et al., 1997, Sheng and Thompson, 2013).Typical applications of PNMs include water retention (Matthews et al., 1995), water permeability (Valvatne and Blunt, 2004), oil wettability (Kovscek et al., 1993, Oren et al., 1998), and capillary pressure (Raoof and Hassanizadeh, 2012).Mechanical coupling in porous materials has been implemented in conjunction with either discrete or finite element models (i.e., recently Yan et al. (2021) and Le et al. (2022)).Such simulations are generally restricted to a relatively narrow range of low degrees of saturation, such as in the range of capillary bridges (Richefeu et al., 2008) and the pendular regime (Scholtes et al., 2009), and only recently were they extended to higher degrees of saturation, such as the funicular regime (Yuan and Chareyre, 2022).Drainage and imbibition paths in clays involving a wetting and a nonwetting fluid induce significant volume changes in clayey or silty geomaterials (Simms andYanful, 2002, Nikooee et al., 2014).As highlighted by Berli and Or (2006), pore deformation upon hydromechanical paths remains understudied, as well as methodologies allowing us to consider any hydromechanical coupling in PNMs (Nikooee et al., 2014).Simms and Yanful (2005) modelled the drying process of compacted clay samples with a deformable pore-network model (DPNM).Their DPNM was designed and validated for a very specific path but has the potential to be extended to more general hydromechanical loading.Simms and Yanful's model is advanced in this paper by i) characterising the deformation of the pores for the case where the fluid pressure in the pores decreases due to the replacement of water under tension with air at atmospheric pressure and ii) identifying a procedure to determine the constitutive parameters (i.e., the stiffness of a single pore upon loading and unloading) based on experimental tests independent to the tests simulated by the DPNM.
The proposed DPNM is based on the binary approach discussed by Pedrotti and Tarantino (2018b), Pedrotti and Tarantino (2018a), Pedrotti and Tarantino (2019), i.e., pores in the clay are assumed to be either fully saturated with water or dry (air-saturated).The DPNM was first tested against the macroscopic and microstructural response of two clays and then used to explain the different mechanisms controlling the volumetric response of clays upon drying.

Materials
Two different clays were chosen for the experimental tests, Speswhite Kaolin (SK), with plastic limit w P = 32 %, liquid limit w L = 64 %, and grain size distribution showing it to have 0.20 silt fraction and 0.80 clay fraction (Fig. 1), and Vitreous China (VC) mixture composed of ball clay (20-28 %), kaolin (25-25 %), quartz (17-31 %), and feldspar (15-23 %).This mixture is currently used in sanitary manufacturing due to its reduced susceptibility to cracking.The plastic limit was found to be w P = 17 %, the liquid limit was w L = 23 %, and the grain size distribution showed it to have a 0.73 silt fraction and a 0.27 clay fraction (Fig. 1).

Drying behaviour
Speswhite kaolin was reconstituted from slurry prepared at a water content equal to 1.5⋅w L and consolidated to 100 kPa vertical stress under saturated conditions in an oedometer cell (80 mm diameter).After consolidation, the sample was unloaded, removed from the cell, and airdried from the top and bottom surfaces to the target water content.Once the target water content was attained, a specimen was cut and trimmed by means of a sampler with a known volume.Suction was measured by means of a high capacity tensiometer (Tarantino and Mongiovi, 2003) at low pore-water tension.A transistor psychrometer was used to measure high pore-water tension (Tarantino and De Col, 2008).Finally, the specimen was placed in an oven at 105 • .As the total volume, V T , the mass of solids, M S , and the water content, w, were measured, the void ratio ( VT − MS/ρ S MS/ρ S ) and degree of saturation ) could be calculated, where ρ W and ρ S are the density of water and the solids, respectively.An additional series of Speswhite kaolin samples were prepared using the same procedure to extract specimens for the pore-size distribution measurements.
Specimens of Vitreous China (VC) were prepared by slip casting in cuboidal plaster moulds (200 mm long, 100 mm wide and 15 mm high), as described in Murray and Tarantino (2018).The suction possessed by the plaster removes water from the liquid slip and consolidates the clay to create a plastic specimen.The slip was prepared at a casting density of 1.84 kg/l, w = 0.35 and added to the plaster mould.The time needed to consolidate the slips of VC with a moisture content of 0.2-0.21 in the mould was approximately 90 min.Once removed from the mould, the bars were air dried to different target water contents.Once the target water content was attained, each bar was sealed and stored for 24 h prior to testing to obtain equilibrium in the water content throughout the bar.From each bar, individual specimens were taken for specific tests (i.e., pore-water tension, void ratio, and pore-size distribution measurements).Pore-water tension was measured by means of a high capacity tensiometer (Tarantino and Mongiovi, 2003) at low pore-water tension.A WP4 dewpoint potentiometer was used at high water tension.

1-D compression test
Compressibility indices for saturated clay and dry powder were derived from oedometer tests on clay reconstituted from slurry (w = 1.5⋅wL ) and dry powder, respectively, for both SK and VC.Clay was consolidated/compressed in steps up to 2220 kPa and then unloaded to 1 kPa.For each step, the saturated clay was allowed to consolidate.The void ratio was back-calculated from the vertical displacement measured in loading/unloading steps and either the final water content (saturated samples) or the final sample volume (dry powder).

Pore-size distribution
The pore size distribution was measured on specimens with water contents ranging from 0.41 to 0.1 for Speswhite kaolin and from 0.18 to 0.10 for Vitreous China.Before mercury intrusion porosimeter (MIP) testing, specimens were dehydrated by means of freeze-drying to preserve their microstructure.Freeze-drying was performed by freezing the specimens using liquid nitrogen and isopentane as described in Tarantino and De Col (2008) and Pedrotti (2016), and the ice water was removed via sublimation in a ScanVac CoolSafe freeze-dryer.
MIP tests were carried out using a Poremaster-60 manufactured by Quantachrome Instruments.The equipment is designed to measure pore entrance diameters in the range from approximately 1000 μm to 0.003 μm corresponding to 1.5-420,000 kPa mercury intrusion pressure.For pressures below 340 kPa, the machine operates on a station different from the one used for high pressures.Blank tests were performed for the two porosimeters to correct errors associated with system compressibility and temperature effects.The intrusion rate was selected to minimise the temperature increase in the system (Simms and Yanful, 2004).

Experimental results
Fig. 2 shows the water retention curves in terms of void ratio, degree of saturation and water content for both Speswhite Kaolin (Fig. 2a and Fig. 2b) and Vitreous China (Fig. 2c and Fig. 2d).
From a mechanical perspective, the volume change presents a relatively complex pattern.As the samples dry out from a slurry state, suction develops, and the void ratio decreases.In the first stage, the clay remains saturated (Fig. 2a and Fig. 2c), and the void ratio decreases linearly with suction (in a semilogarithmic space).
Desaturation (Air Entry Value -AEV) occurs at a pore-water tension of approximately 1000 kPa and water content of approximately 0.32 for SK (AEV in Fig. 2b) and a pore-water tension of approximately 800 kPa and water content of approximately 0.13 for VC (AEV in Fig. 2d).As the clay desaturates, the void ratio decreases at a very low rate, reaches a minimum, and then increases as the pore-water tension increases.
In terms of hydraulic behaviour, the clay moves from the saturated state to the funicular and eventually the pendular state, i.e., almost the entire range of degrees of saturation is explored upon the drying process.Fig. 3 shows the evolution of the cumulative pore size (Fig. 3a and Fig. 3b) and frequency pore size distribution (Fig. 3c and Fig. 3d) upon drying for SK and VC, respectively.The pore-size distribution was measured on 6 samples for SK and 4 samples for VC at different water contents (solid circles in Fig. 2b and Fig. 2d).
For SK (Fig. 3a), the cumulative void ratio decreases from 1.26 to 0.75 as the sample dries out, and the water content decreases from 0.18 to 0.13.As the water content further decreases from 0.28 to 0.09, the change in the cumulative void ratio is reversed (the cumulative void ratio attains a value of approximately 0.85 at the end of the drying process).Similarly, the cumulative void ratio for VC decreases from 0.39 to 0.31 as the water content decreases from 0.18 to 0.13 (Fig. 3b).As the water content further decreases to 0.10, the change in the cumulative void ratio is reversed (the cumulative void ratio attains a value of approximately 0.32 at the end of the drying process), although not as much as for the SK samples.
For SK, the modal pore size decreases as the samples dry as long as the sample remains saturated (water content from 0.47 to 0.32).Similarly, the void ratio frequency associated with pore sizes larger than 0.3 μm decreases.On the other hand, the modal pore size does not change for samples with water contents smaller than 0.28.However, the void ratio frequency of larger pores tends to increase (see void ratio frequencies in the pore range 0.3-1 μm for samples having water contents equal to 0.09, 0.21 and 0.28, respectively, in Fig. 3c).
A similar trend can also be observed for VC.For water contents from 0.18 to 0.13 (when the sample is saturated), the modal pore size decreases as the samples are dried.On the other hand, the modal value increases as the clay desaturates, as does the frequency of larger pores (Fig. 3d).

A binary approach for microscale modelling of the drying of saturated clays
The experimental characterisation of the microstructure of partially saturated kaolin presented in Pedrotti and Tarantino (2019) showed that the pore size distribution (PSD) of kaolin reconstituted from slurry (aka water-saturated sample in Fig. 4) and the PSD of compressed kaolin dry powder (aka air-saturated sample in Fig. 4) overlap with the modal size associated with the micro-and macropores of the compacted samples shown in Fig. 4. Based on the assumptions that compacted clay is made of a water-saturated part (pores filled with water) and an air-saturated part (pores filled with air), the following hypothesis was proposed for partially saturated states: -Water-saturated pores behave like water-saturated clay -Air-saturated pores behave like dry clay The above hypothesis implies that two different effective stresses are defined for the water-saturated and air-saturated subregions, and accordingly, two difference sets of constitutive relationships are needed to model the deformation of the water-saturated and air-saturated subregions.
A conceptual model based on the above hypotheses was formulated and successfully used to simulate the macroscopic behaviour of hydromechanical paths on partially saturated kaolin, i.e., compaction, wetting + reloading, drying, and direct shear (Pedrotti andTarantino, 2018a, Pedrotti andTarantino, 2019).
This binary approach is extended here to the pore scale: -Pores that are water-saturated reduce their volume as suction increases.In the saturated range (A to B in Fig. 5), all pores reduce their size -Pores that desaturate are subjected to a rebound because water under tension is replaced by air at atmospheric pressure.In the unsaturated range (steps B to C in Fig. 5), some pores desaturate and, hence, rebound.The pores that remain saturated will continue to decrease in size as pore-water tension increases.
According to the hypotheses presented above, when the pore is water-saturated and suction increases, the effective stress for the given pore will be calculated as σ i − u i , where σ i is the total stress acting on the given pore i, and u i is the pore pressure within the pore.The compressibility of the pore will be macroscopically associated with the compression indices λ w or k w , i.e., the slopes of the normal consolidation line (NCL) and unloading-reloading line (URL) for the clay reconstituted from slurry, respectively.When pore i desaturates, the effective stress will be calculated as σ i (the air pressure is assumed to be zero).The compressibility of the airsaturated pore will be macroscopically associated with either λ a or k a , which are the compression indices of the NCL and the URL for the dry clay, respectively.
If water is present within a pore, the pore shrinks; if the pore desaturates, a rebound mechanism is triggered owing to the decrease in effective stress acting on the pore.These competing mechanisms (pore shrinkage and pore rebound) might result in either a decrease or an increase in the void ratio at the macroscopic level depending on the number and size of pores that shrink or rebound.

|Pore geometry
Pore network models for deformable clayey materials have been little explored thus far (Simms and Yanful, 2002, Simms and Yanful, 2005, Rostami et al., 2013, Khaksar et al., 2013).A number of challenges need to be addressed, including hydromechanical coupling at the pore scale and the characterisation of the pore space.Pore network models have been generally designed for sand-like porous materials (having particles of quasispherical shape), where the pore space can generally be described by a system of pores and channels.Furthermore, the pore geometry and connectivity in sand-like porous materials can be directly visualised by microscopic imaging methods, e.g., X-ray computed tomography (amongst many Blunt et al. (2013)).
On the other hand, clay particles are generally plate-like, and their surface is electrically charged.This generates complex particle configurations depending on the electromechanical interaction between clay particles and pore water (Pedrotti andTarantino, 2018b, Pedrotti andTarantino, 2019).
Clay particle sizes fall in the range of microns, thus creating pores in the submicron range.This makes it virtually impossible to directly visualise the 3-D pore space by means of current microscopic techniques.Due to the lack of direct experimental observations, the 3-D geometry of the pore space can only be speculated.As discussed in the Appendix A, reasonable and convenient pore geometry can be represented by pores that are directly connected to each other without interconnecting channels.
The aim of the DPNM is to translate the hydromechanical interactions at the pore scale inferred from the pore size distributions and to test its capability to quantitatively reproduce the macroscopic behaviour observed in terms of void ratio (volume change) and water content (water mass change).
The pore space was defined as the volume existing between two parallel plates at distance D i and surface area A (Fig. 6).This was considered appropriate for a clayey geomaterial and is also in line with the approach presented by Simms and Yanful (2005).The void volume, V i V , for each pore can be expressed as follows: The network model idealises the porous medium as a threedimensional lattice of intersecting pores of different sizes.Each pore is connected with 26 neighbouring pores (Fig. 7) in 3-D space.

Pore hydraulic behaviour
The process of pore desaturation for the pores at the interface between air and water is assumed to be controlled by the Young-Laplace equation.The difference between the pressure of the water (wetting fluid), P w,i , and the pressure of the air (nonwetting fluid), P a , for a meniscus having a surface with radii of curvature R 1,i and R 2,i, respectively, is given by: where γ is the air-water surface tension.By assuming that the single pore is characterised by two parallel plates and that the plate distance D i is small compared to the plate dimension, the radius R 2,i becomes very large.The air-water interface turns into a cylindrical surface characterised by a radius of curvature R 1,i defined by: where θ is the contact angle.The Young-Laplace equation can be written as a function of D i as follows: The fluid state of each pore is defined in a binary fashion, i.e., the pore can be either water-saturated or dry (air-saturated).No partially saturated pores are considered.
As the clay geomaterial is dried from one of its outer boundaries and the pore-water pressure decreases (i.e., becomes more negative), a saturated pore becomes dry if the following criteria are met: -The water-saturated pore is connected to the external air via a continuous path of dry pores.-The pore-water pressure becomes lower than the desaturation pressure given by Equation (4).

Pore mechanical behaviour
Following the approach adopted by Simms and Yanful (2005), the pore-scale strain ε i of pore i is assumed to be proportional to the change in pore-scale effective stress Δp′: where Δd i is the change in the size of pore i, d i is the pore dimension, and C is a stiffness parameter with dimensions of 1/Pa.Following such an assumption, it is shown later in the paper that the strain is uniform for all the pores and is not dependent on the pore size itself.
The strain at the pore scale was hence assumed to be uniform for all the pores (as per Equation ( 5)) and equal to the macroscopic strain ε, as shown schematically in Fig. 8: This is the simplest assumption about the distribution of strain across the pores and was favoured with respect to any other arbitrary assumption about strain distribution.It is shown later in the paper that such an assumption was adequate to simulate the evolution of pore-size distribution and volume upon drying.
The strain at the pore scale was assumed to be controlled by the effective stress acting on the pore, σ ′ i , in turn given by the difference between the total stress acting on the pore, σ i , and the pressure of the fluid in the pore, P i .
The pore total stress σ i (as represented in Fig. 8) was assumed to be equal to the external isotropic total stress σ, i.e., the external stress was assumed to diffuse uniformly across the pores.Again, this simple assumption was favoured with respect to any other arbitrary assumption about the distribution of stresses across the pores.The fluid pressure, P i , in the pore was assumed to be equal to either the external air pressure P A (air-saturated pore) or the pore-water pressure P W (water-saturated pore).The pore-air pressure and pore-water pressure were assumed to be in equilibrium (i.e., uniformly distributed) across the air-saturated and water-saturated pores, respectively.In summary: To derive the relationship between stress and strain at the pore scale, the relationship between stress and strain at the macroscale was first considered and then scaled down to the pore scale, as illustrated below.

Macroscopic constitutive law
For the case of saturated clay, the volume change behaviour is controlled by the macroscale effective stress σ′, the difference between the total stress σ and the pore-water pressure u w (Terzaghi, 1936) Volume change typically shows elastic and plastic components.Fig. 9 shows the evolution of the void ratio e (void volume to solids volume ratio) with the effective stress upon loading and subsequent unloading.The void ratio decreases when the clay is loaded from a virgin state.The relationship between the void ratio and effective stress is generally described by a semilogarithmic function where λ W is the compressibility upon virgin loading.
If the effective stress is reversed, only a part of the deformation accumulated upon virgin loading is recovered, i.e., plastic deformation occurs.Again, the relationship between the void ratio and effective stress upon an unloading path is generally adequately described by a semilogarithmic function: where κ W is the compressibility upon unloading-reloading (κ W < λ W Eq. ( 11) applies when the current effective stresses are lower than the maximum ever experienced (overconsolidated state).
The equations used to represent the compressibility of clay (Eq.(10) and Eq. ( 11)) cannot be extrapolated to zero effective stress (void ratio goes to infinity).The assumption made in this paper is that the void ratio attains a constant value at a threshold stress σ 0 ′ (Fig. 9).
For the case of clay powder (air-saturated clay), Pedrotti and Tarantino (2018a) and Pedrotti and Tarantino (2019) showed that the same constitutive laws also apply provided that different constitutive parameters are used, λ A and κ A for the virgin and unloading-reloading compression lines, respectively (Fig. 9):

Pore-scale constitutive law
Clay prepared from slurry is in a 'virgin' state.If subjected to drying, the pore-water pressure decreases (i.e., becomes more negative), the effective stress increases, and the clay moves along the virgin compression line characterised by the virgin compressibility λ W (Fig. 9).This process is also assumed to occur at the pore scale, and the mechanical constitutive law for a saturated 'virgin' pore subjected to negative pore water pressures was therefore derived from the virgin compression line.
Let us first consider the volume of the voids, V V .This is the cumulative sum of the volume of voids associated with each pore, V Vi , and according to Equation (1), we can write: By introducing the average pore size D as: and differentiating Equation ( 14), we can write On the other hand, let us rewrite Equation ( 10) by making explicit the change in the volume of voids, dV V , and the volume of the solids, V s : By combining Equations ( 16) and ( 17) and by dividing the right-and left-hand sides by the average pore size D, we obtain where ε is the macroscopic volumetric strain.
Under these assumptions, the pore volumetric strain can be written as The number of pores N, the volume of voids V S , and the particle area A are not known a priori.However, they can be derived from the cumulative pore-size distribution (PSD), as shown in Fig. 10.Let us divide the pore-size range into several classes, N class .The total number of pores, N, is then given by where n j is the number of pores for the given class j.The volume of voids associated with each pore class, V V j , can be written as where D j is the size characterising the given class j.The term A N VS can then be derived using Equations ( 20) and ( 21) as follows: where Δe j is the void ratio associated with each pore class j, as shown in Fig. 10.
Finally, the pore strain ε i can be rewritten as a function of the increment of pore effective stress as follows.
The term in the square bracket represents the stiffness of the watersaturated pore upon virgin loading.In this term, the only parameter that needs to be calibrated is the water-saturated compressibility index, λ W , which can be derived from the macroscopic compressibility of the watersaturated sample.
When an initially saturated pore desaturates (i.e., it becomes filled with air), its effective stress decreases, and the pore rebounds.The pore lies on an unloading-reloading line associated with dry (air-saturated) pores, and the mechanical constitutive law for a dry pore subjected to negative pore-water pressures can therefore be derived from the unloading-reloading line of air-saturated clay.Using an approach similar to the one adopted to derive equation ( 23), the following constitutive law can be derived for the pores that desaturate: Again, the term in the square bracket represents the stiffness of the pore upon unloading-reloading for an air-saturated sample.Again, the only parameter that needs to be calibrated is the loading-reloading airsaturated compressibility, κ A , which can be derived from the macroscopic compressibility of the air-saturated sample.

Initial conditions and input parameters
The pore network model was implemented from scratch in Python and is available at the following repository (https://github.com/PedroMat8/DPNM_Computers-and-Geotechnics).The initial distribution of the pores in the model was generated from the PSD of the sample dried to a water content of 0.47 for SK and 0.18 for VC.Fig. 11 shows the experimental PSD and the initial PSD generated for the DPNM.
For the simulation, a 3-D network of 125,000 pores (50 × 50 × 50) with a coordination number of 26 was adopted.In the model, the water surface tension, γ (), was assumed to be equal to 0.072 N/m, and the receding contact angle, θ ( • ), was set equal to 0. Pores were assumed to be initially fully saturated with water.The clay particle specific gravity was assumed to be 2.605 for both clays.The effect of pore uniformity on the PSD was not considered, as the investigated hydromechanical path was a drying process, hence equivalent to the intrusion of nonwetting mercury into dry soil (Romero and Simms, 2008).

Calibration
Compressibility indices were calculated from a 1-D consolidation performed in an oedometer cell on reconstituted and dry powder samples for both SK and VC clays.Fig. 12 shows the experimental data of the 1-D compression test (dots) and the fitting curves (solid lines) used to interpolate the loading and unloading curve segment.The two compression indices, λ W and κ A, are also reported in the figure for both SK and VC.The compression indices λ W and κ A derived from the 1D compression test were assumed to coincide with the compression indices under isotropic conditions, i.e., the earth coefficient at rest k 0 was assumed to be constant.It is acknowledged that this assumption may not hold for overconsolidated dry powder; however, the error introduced is considered to be minor in light of the correction proposed in Fig. 9.

Sensitivity analysis for the number of pores and coordination number
A sensitivity analysis was carried out to validate the number of pores and coordination number adopted in the analyses.Sensitivity analysis was performed for the case of SK.
The air-entry point in terms of water content and matric suction was considered to benchmark the results from the simulations (Fig. 13a).Fig. 13b shows the effect of the cubic matrix dimension on the water content at the air-entry value.Cubic matrices ranging from 10 pores per side (10 3 pores in total) to 100 pores per side (10 6 pores in total) were initially investigated (using a 3-D configuration with a coordination number equal to 26).The simulated water contents at the air entry value (w AEV ) are shown in Fig. 13b, where they are compared with the experimental value.It can be observed that beyond a matrix size larger than 50 pores per side, the simulated w AEV fluctuates around the experimental value.
To investigate the effect of the coordination number and the pore space distribution (2-D or 3-D configuration), the results of four different analyses were compared.Two analyses were carried out on a square matrix (2D) with 300 pores per side (90,000 pores in total) but coordination numbers of 4 and 8 (open symbols in Fig. 13c).Similarly, two analyses were carried out on a cubic matrix (3-D) with 50 pores per side (125,000 pores in total) but coordination numbers of 6 and 26 (solid symbols in Fig. 13c).The simulated values of suction at the entry value are compared with the experimental value in Fig. 13c.Even if no simulations returned exactly the experimental value, the higher the coordination number, the smaller the error for both the cases of 2D and 3D spatial distribution.The high coordination number needed is expected if one considers the open structure that characterises clayey material, where the volume of voids is often higher than the volume of solids (i.e., void ratio higher than 1).

Stress path and boundary conditions
Water contents of 0.47 for SK and 0.18 for VC were considered the initial conditions for the drying test.Initial pore-water tensions were set equal to 99 kPa for SK and 104 kPa for VC, and no external total stress was applied.Under these conditions, the sample was fully saturated.To simulate the drying process, the pore-water tension was increased in incremental steps until all pores became dry (air-saturated).Pores at the 6 outer boundaries were assumed to be directly in contact with air and hence free to desaturate if the imposed suction was higher than the suction needed to desaturate the pore.Upon each step of pore-water tension increase, the hydraulic status of the pore (saturated or dry) was updated according to Equation (4).For each step, the pore system was iteratively scanned until no new pores were desaturated and the system could be considered to be at equilibrium.
Pores that remain saturated are subjected to an increase in effective stress since pore-water tension increases.Their volume therefore decreases according to Equation ( 23).Pores that desaturate are subjected to a decrease in effective stress down to zero (associated with the zero pressure caused by the air filling the pore).Their volume therefore increases according to Equation (24).To calculate the increase in volume of the pores that desaturate, the threshold effective stress σ 0 ′ was set equal to 1 kPa.
Fig. 14 shows a representation of the pore matrix at different moisture contents for the case of SK.To ease the visualisation, a 50 × 50 × 50 pore network is shown.Matrix cells corresponding to water-saturated pores are represented (blue in Fig. 14), and dry pores are transparent.As expected, the pores desaturate as the pore-water tension increases.

Simulation of volume change and water retention behaviour
The PSD at the onset of the drying process measured by the MIP was used to set the initial condition for the DPNM.As a result, the initial void ratio simulated by the DPNM coincides with the cumulative pore volume intruded by mercury.
There is generally a small discrepancy between the volume of the pores in the clay sample determined experimentally via the measurement of the total volume of the specimen and the volume intruded by the mercury (as the mercury does not intrude the very smallest pores).For this test, the discrepancy in void ratio resulted in Δe = 0.098 for SK and Δe = 0.102 for VC.This difference was used to correct the void ratio derived from the DPNM simulation.
Fig. 15 compares the experimental results with the simulation.Fig. 15a and Fig. 15c show the change in void ratio against the porewater tension for SK and VC, respectively, whereas Fig. 15b and Fig. 15d show the change in water content against pore-water tension, respectively.The experimental data are simulated satisfactorily in terms of both void ratio and water content.The simulation of the water retention curve deviates slightly from the experimental data at very high pore-water tension for the case of SK.Such a difference might be due to mechanisms of water retention different from capillarity.
Notably, the complex volume change behaviour occurring upon desaturation (initial shrinkage followed by a rebound) is well simulated by the DPNM.Upon desaturation, the volume of the pores that remain saturated continues to decrease in size, whereas the pores that dry out rebound.The overall volume change is the result of the competition between these two processes.At the onset of desaturation, most of the pores are saturated, and shrinking prevails.As the clay dries out, more and more pores desaturate, and an overall rebound is observed.

Prediction of the evolution of PSD
Fig. 16 and Fig. 17 compare the PSD obtained from the simulation with the experimental PSD for the tests reported in Fig. 3 for SK and VC, respectively.The model shows a satisfactory prediction for water contents spanning from 0.47 to 0.09 for SK and from 0.18 to 0.10 for VC.This further confirms the capability of the DPNM to capture the complex coupled hydromechanical behaviour of clay upon drying.

Conclusions
The paper has presented a simple pore-network model to elucidate the fundamental pore-scale mechanisms controlling the hydraulic and volumetric response of clay upon air-drying.
Two clay materials have been tested, Vitreous China and Speswhite Kaolin.The hydromechanical response of the clay upon air-drying drying was characterised in terms of both moisture content (mass of water per mass of solids) and void ratio (volume of voids per volume of solids) versus pore-water tension (capillary suction).The evolution of the void ratio upon drying appears to be characterised by three stages.As porewater tension increases, the void ratio initially decreases considerably while the clay remains in a saturated state.The void ratio then tends to level off as soon as the clay desaturates and eventually increases again at relatively high pore-water tension.The macroscopic measurements were complemented by the measurement of the pore-size distribution via mercury intrusion porosimetry at the different points along the drying curve.The change in void ratio accumulated in the saturated range appears to be associated with a significant reduction in the larger pores.In contrast, when the clay desaturates, some of these larger pores clearly appear to rebound.
A conceptual micromechanical model was then developed to explain different stages in the volumetric response of clay upon air-drying.In the saturated state, all the pores shrink as pore-water tension increases.When the clay desaturates, the larger pores dry out and rebound, whereas the smaller pores remain saturated and still shrink.These two competing mechanisms are responsible for the volumetric response of the clay following desaturation.
Based on this conceptual micromechanical model, a deformable pore-network model (DPNM) was formulated to allow quantitative coupling of the hydraulic and mechanical responses of clay at the pore scale upon unconstrained drying.This model was purposely maintained very simply and only requires the initial pore-size distribution and two independently determined constitutive parameters as inputs.The porenetwork model was successfully validated against experimental data at the micro-and macroscale, thus corroborating the assumptions about the competing pore-scale mechanisms controlling the mechanical response of clay upon air-drying.The presented model has the advantage of requiring reduced input parameters (one initial pore size distribution and compressibility of clay in dry and reconstituted states).If properly modified by introducing strain compatibility conditions, the PNM can potentially be used to simulate cracking.

Credit author statement
Matteo Pedrotti is the corresponding author and wrote the paper.He contributed to the paper with the conceptualization of the idea, the experimental work on kaolinite, developed the numerical code and run the simulations.
Long Xu run the initial preliminary numerical simulations.Ian Murray carried out the experimental work in Vitreous Clay.
Alessandro Tarantino supervised the whole work, contributed to the conceptualization of the paper and in writing the final version.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
imaging methods are currently available for wet clays and that there is therefore no experimental evidence of pore-space configuration in clays.

Fig. 2 .
Fig. 2. Drying behaviour macroscopic behaviour.a) Void ratio and degree of saturation vs. pore-water tension for Speswhite kaolin.b) Water content vs. pore-water tension for Speswhite kaolin, c) void ratio and degree of saturation vs. pore-water tension for vitreous China and d) water content vs. pore-water tension for vitreous China.Samples tested for MIP are reported as solid circles.

Fig. 10 .
Fig. 10.Derivation of the PSD characteristic A N VS .

Fig. 12 .
Fig. 12. Compressibility indices derived from the 1-D compression test.a) SK reconstituted from slurry, b) SK powder, c) VC reconstituted from slurry and d) VC powder.

Fig. 13 .
Fig. 13.Sensitivity analysis.a) Air entry value, b) Effect of cubic matrix dimension and c) Effect of coordination number.

Fig. 14 .
Fig. 14.DPNM central slice.Top left: Initial pore size.From top left to bottom right: Saturation of the pores forming the DPNM for moisture contents ranging from 0.4 to 0.01.

Fig. 15 .
Fig. 15.Simulation of volume change and water retention behaviour.