Network modelling of fluid retention behaviour in unsaturated soils

The paper describes discrete modelling of the retention behaviour of unsaturated porous materials. A network approach is used within a statistical volume element (SVE), suitable for subsequent use in hydro-mechanical analysis and incorporation within multi-scale numerical modelling. The soil pore structure is modelled by a network of cylindrical pipes connecting spheres, with the spheres representing soil voids and the pipes representing interconnecting throats. The locations of pipes and spheres are determined by a Voronoi tessellation of the domain. Original aspects of the modelling include a form of periodic boundary condition implementation applied for the first time to this type of network, a new pore volume scaling technique to provide more realistic modelling and a new procedure for initiating drying or wetting paths in a network model employing periodic boundary conditions. Model simulations, employing two linear cumulative probability distributions to represent the distributions of sphere and pipe radii, are presented for the retention behaviour reported from a mercury porosimetry test on a sandstone.


Introduction
Realistic modelling of fluid retention behaviour of porous materials is important in many fields.Selection of parameter values within conventional analytical expressions for retention curves, such as [1], is normally based on back-fitting experimental retention curves.This does not provide insight to improve the understanding of the influence of microstructure on the fluid storage properties of the material.Discrete numerical modelling of fluid storage and transport in porous materials at the microstructural scale can improve this understanding.
Since the pioneering work of [2], network models have been widely used for fluid retention and hydraulic conductivity modelling of porous materials, particularly in fields other than geotechnical engineering.Authors such as [3] and [4] investigated the influence of mean coordination number, sphere and pipe size distributions and spatial correlation on retention curves and relative permeability curves, including the influence of phenomena such as trapping and snap off.
[2]- [4] aimed to improve qualitative understanding of the influence of microstructure on macroscopic fluid retention and transport properties of porous materials.Subsequently [5] proposed a calibration technique to quantitatively match experimental mercury intrusion and extrusion curves.This approach managed to reproduce retention curves for values of degree of saturation between 0.25 and 0.6.This calibration approach was then extended by [6], where the main objective was to deduce the material microstructure from experimental fluid retention curves and then evaluate the conductivity curves.Later this model was used in [7] for investigation of CO 2 storage capability of sandstones due to trapping.
Previous studies that attempted to quantitatively reproduce retention properties of porous materials by network modelling used computationally demanding algorithms and also resulted in relatively poor matching towards the tails of retention curves.Also, most models involved the generation of 3-D networks based on cubic arrays that induce a structural regularity that could produce problems for potential usage in hydromechanical analyses.Also periodic boundary conditions (PBCs) were applied on only 4 lateral faces of a statistical volume element (SVE), so that fluid could be injected from one of the remaining faces.
In the present work, a novel computationally efficient irregular network model for achieving quantitatively realistic fluid retention modelling is described.Original features of the model include a new way of setting up analyses with full PBCs and a novel pore volume scaling technique, which avoids the need for excessively large networks to represent materials with a wide range of void sizes.The approach is suitable for future extension to hydro-mechanical and multi-scale numerical modelling.

Network model
The discrete numerical model for the analysis of fluid retention is based on a three-dimensional network model of one-dimensional pipes that connect spheres, with the spheres representing soil voids and the pipes representing inter-connecting throats.
The network construction is based on a dual Delaunay-Voronoi tessellation.The domain is initially saturated with points placed randomly but with a restriction on the minimum separation distance d min .The E 2016 -Delaunay tessellation (consisting of irregular tetrahedra with their vertices at the points) is then applied and consequently its dual Voronoi tessellation is generated (consisting of irregular polyhedra, with each face of a polyhedron perpendicular to a Delaunay edge and intersecting this Delaunay edge at its mid-point).The spheres are then placed on the polyhedra vertices and the pipes along the polyhedra edges.This is a suitable approach for subsequent hydro-mechanical analyses, where the polyhedra are idealisations of the porous material grains surrounded by voids and throats, as discussed in [8].Grain interactions can be modelled by placing mechanical elements along the Delaunay edges.A schematic presentation of four Voronoi polyhedra and their corresponding Delaunay edges is presented in figure 1, with each polyhedral grain scaled down in the figure in order to provide a clearer visual impression (in the actual network the grains are in contact).The work presented in this paper was restricted to analysis of fluid retention behaviour in a non-deforming material i.e. it did not include mechanical elements.
Spheres and pipes each have a distribution of radii, each with a specified probability distribution function.Sphere radii are allocated first, on a random basis, with no spatial correlation.Pipe radii are then allocated to locations, with the smallest radii pipes allocated first, commencing with locations connecting to the smallest spheres.

Periodic boundary conditions
The numerical analysis of fluid retention presented in this paper was conducted on a network representing a cubical statistical volume element (SVE).Periodic boundary conditions (PBCs) were applied on all 6 faces of the SVE, so that the SVE was considered to repeat indefinitely in all directions.Traditional techniques for PBCs, of Lagrange multipliers or four driving nodes, are not well suited for irregular Voronoi and Delaunay networks.Instead a different PBC technique was used.This technique was presented in [9] for mechanical lattice elements in a 2-D domain, and in the current work it was extended to 3-D transport networks.This PBC approach requires special treatment of the network construction, because the elements are allowed to cross the SVE boundaries.An example of a 3-D periodic transport network is presented in figure 2. It can be observed that the pipes, presented in grey, cross the SVE boundaries indicated by the black edges of the cell.

Pore volume scaling
For a material with a wide range of void sizes, a SVE would require a very large number of smaller voids in order to also include a statistically representative number of larger voids.With a standard network, this would require a SVE with an impractically large number of spheres.In addition, a standard network could not capture the fact that the average centre-to-centre spacing of smaller voids should be much less than the corresponding spacing for larger voids, otherwise adjacent small voids would be unreasonably far apart (suggesting an exceptionally low local porosity) and two adjacent large voids would be impossibly close together (suggesting overlap of the two voids).In order to avoid these problems, a new approach was introduced, where each sphere within the network represents a scaled number of voids, with this scaling number increasing as the sphere size reduces.
In the pore volume scaling technique, each sphere of radius r s is taken to represent not a single void but a non-integer number N of voids, each of radius r s , where N is given by where r sm is the mean radius of the spheres.Eq.1 means that each sphere of radius smaller than the mean represents more than a single void, whereas each sphere of radius greater than the mean represents less than one void.This scaled number of voids means that each sphere represents a volume V v of voids given by Eq.2 shows that each sphere, whatever its radius, represents the same volume of voids.The volume of pipes is ignored in calculating the porosity or degree of The advantage of the pore volume scaling approach is that the network model is able to represent a SVE containing a large number of small voids without excessive computational effort.In addition, local values of porosity internally within the SVE are realistic, rather than being unreasonably low in regions around smaller spheres and impossibly high in regions around the largest spheres.A disadvantage of the pore volume scaling approach is, however, that modelling of connectivity between voids, which is always imperfectly represented in a network model, may be represented less realistically than without scaling.This can be illustrated by considering the case of a sphere that is significantly smaller than the mean sphere radius r sm..This sphere represents a substantial number N of voids, rather than a single void.The scaling technique means that local connectivity between these N voids is overstated, as they are all assumed to be at the same location in the network and are therefore perfectly connected to each other.In contrast, external connectivity between these N voids and other voids is understated, as only 4 pipes connect the single sphere representing all these N voids to other spheres.Reverse arguments apply to voids that are larger than the mean radius.

Modelling retention behaviour
Capillary suction P c , applied uniformly throughout the SVE, was the driving variable for drying and wetting of the SVE during modelling of retention behaviour.The applied value of P c was monotonically increased during a drying process and monotonically decreased during a wetting process.

Representing drying and wetting processes
The rules that determine, for a given value of P c , which spheres are filled with the wetting fluid and which are filled with the drying fluid are based on the Young-Laplace equation where γ is the surface tension and θ is the contact angle formed between the fluid-fluid interface and the solid (measured on the wetting fluid side).r is the radius of a during a drying process and the radius of a sphere during a wetting process, to represent the so-called "ink bottle" effect.A drying or wetting process is modelled by considering the intruding fluid moving from a void already filled with that fluid to a connected void not previously filled with that fluid i.e. direct connection to a void already filled with the intruding fluid is imposed as a requirement for a void to fill with the intruding fluid.The drying fluid is the intruding fluid during a drying path, whereas the wetting fluid is the intruding fluid during a wetting path.
Figure 3 illustrates the situation during a drying process.The hatched area corresponds to the wetting fluid.The drying fluid has already intruded through pipe 1 into sphere A at a value of P c corresponding to the application of Eq.4 with r as the radius of pipe 1.As the value of P c is increased further during the drying process, the next consideration is when the drying fluid will intrude further from sphere A along either pipe 2 or pipe 3. Pipe 2 is of larger radius than pipe 3, and therefore the first action to occur is intrusion of the drying fluid along pipe 2 and into sphere B, at a value of P c corresponding to the application of Eq.4 with r as the radius of pipe 2. If either pipe 4 or pipe 6 was larger than pipe 2 the intrusion would continue along this pipe into an additional sphere without need for further increase of P c . Figure 4 illustrates an equivalent situation for a wetting process, with the hatched area again representing the wetting fluid.Wetting fluid has already entered sphere A from pipe 1, at a value of P c corresponding to the application of Eq.4 with r as the radius of sphere A, and immediately moved into pipes 2 and 3.As the value of P c is reduced further during the wetting process, the next consideration is when the wetting fluid will intrude into sphere B or sphere C. As sphere C is smaller than sphere B, the first thing that happens is intrusion of the wetting fluid into sphere C, at a value of P c corresponding to the application of Eq.4 with r as the radius of sphere C, and immediate further movement of the wetting fluid into pipes 5 and 7.If either pipe 5 or pipe 7 connected to an additional sphere that was smaller than sphere C, the wetting fluid would continue into this sphere without need for further decrease of P c .
At each value of P c during a drying or wetting process, the degree saturation S r (expressed in terms of the wetting fluid) is calculated by considering the proportion of spheres filled with the wetting fluid (bearing in mind that each sphere, whatever its radius, represents the same volume of voids).Hence, if m is the total number of spheres in the SVE and m w is the number of these spheres filled with the wetting fluid, the degree of saturation In analyzing drying and wetting processes, no possibility of trapping of the non-intruding fluid is considered, because there is no consideration of whether there is appropriate connectivity to provide an exit route for the non-intruding fluid (connectivity requirements are applied to the intruding fluid, but not to the non-intruding fluid).In practice, trapping of the non-intruding fluid is not an issue if the non-intruding fluid can escape by diffusing through the intruding fluid and provided that drying or wetting is performed sufficiently slowly for the diffusion process to dissipate excess pressure in the trapped non-intruding fluid.

Initiation
With the use of PBCs, it is not possible to start a drying path from a fully saturated condition (S r = 1) or a wetting path from a fully dry condition (S r = 0), because there is no external boundary from which to introduce the intruding fluid.Instead, a drying path must start with at least one sphere within the SVE already filled with the drying fluid and a wetting path must start with at least one sphere already filled with the wetting fluid.
Investigation showed that realistic modelling of a drying or wetting path could not be achieved by starting with only a single sphere filled with the intruding fluid, because of the extremely low connectivity (only 4 pipes) from a single sphere and the consequent possibility that extreme changes of P c might be required for the intruding fluid to break out from this first sphere (dependent on only the radii of the 4 pipes in the case of a drying path or the radii of the 4 spheres on the other ends of these pipes in the case of a wetting path) or to break out from the local region immediately surrounding the first sphere.In order to achieve realistic modelling, drying or wetting paths should start with a small number of spheres already filled with the intruding fluid and these "seeding" spheres should not simply be selected at random.Instead, to achieve realistic modelling of a drying path, it should always be preceded by modelling of a previous wetting path finishing with a small number of spheres (the seeding spheres) still filled with the drying fluid.A similar condition applies for realistic modelling of a wetting path, which should always be preceded by modelling of a previous drying path to leave an appropriate number of seeding spheres filled with the wetting fluid.
Investigation of this initiation process allowed an appropriate value to be selected for the number of seeding spheres (and hence the appropriate maximum starting value of S r for realistic modelling of a drying path or appropriate minimum starting value of S r for realistic modelling of a wetting path).This was confirmed by the fact that on subsequent cycling between these maximum and minimum values of S r each cycle concluded with exactly the same configuration of spheres filled with the drying fluid.

Application
The network model, including the pore volume scaling technique and the initiation procedure, both described above, was used to study the retention behaviour reported by [5] from a mercury porosimetry test on sandstone.

Procedure and input parameters
In a mercury porosimetry test, air is the wetting fluid and mercury is the drying fluid, with mercury intrusion and extrusion corresponding to drying and wetting paths respectively.The contact angle θ measured on the air (wetting fluid) side is 40 o (see Table 1).
The square data points in figure 5 show the experimental results of the wetting (mercury extrusion) and drying (mercury intrusion) curves respectively in S r -P c space, where S r is the degree of saturation of the air (the wetting fluid).For the wetting path, the maximum value of S r reported in the experimental results was 0.6, because the authors in [5] stated that trapping of mercury prevented significant further increase of S r , presumably because the wetting (mercury extrusion) was performed insufficiently slowly to allow escape of trapped mercury by evaporation and subsequent diffusion through the air.
With the pore volume scaling technique described in Section 2.2, the fact that each sphere represents the same volume of voids means that the drying curve predicted by the network model (in S r -P c space) depends upon only the distribution of pipe (throat) sizes and the network connectivity i.e. it does not depend on the distribution of sphere (void) sizes.Similarly, the wetting curve predicted by the network model depends upon only the distribution of sphere sizes and the network connectivity i.e. it does not depend upon the distribution of pipe sizes.Probability distribution functions for the pipe and sphere radii were therefore selected by interpretation of the experimental results of the drying and wetting curves respectively.
Figure 6 shows the experimental results re-plotted as S r against back-calculated sphere or pipe radius r, by using Eq.4 to convert values of P c to values of r.Of course, even for the idealised case of a network model simulation where the pore volume scaling technique was employed, these plots of S r against back-calculated values of r would not represent the true cumulative probability distribution of sphere or pipe radii, because of the influence of network connectivity.For example, a network model wetting curve simulation presented in this way would only give the cumulative probability distribution of sphere sizes used in the network if there had been perfect connectivity, with access of the wetting fluid to all spheres throughout the wetting process, so that spheres could be filled with the wetting fluid in the exact sequence of increasing sphere size.Similarly, a network model drying curve simulation of S r against backcalculated value of r would only correspond to the true cumulative probability distribution of pipe radii used in the network if network connectivity allowed filling of pipes in the exact sequence of decreasing pipe size and if each pipe then provided access to a sphere that was not previously filled with the drying fluid (the latter is not true because some spheres will already have been filled with drying fluid entering by another route).Despite these issues, it is informative to use the experimental results of S r against back-calculated radius, shown in figure 6, to select the distributions of sphere and pipe radii for the network model.6 shows that the experimental plots of S r against back-calculated value of r (from Eq.4) for both drying and wetting paths can each be approximated by a straight line, at least up to a certain value of S r .For the drying curve this is true up to approximately S r = 0.8 and for the wetting curve it is true up to the final experimental point (S r = 0.6).Two best-fit straight lines were fitted to these parts of the two experimental sets of data, as shown in figure 6, to give an equation for the cumulative probability distribution of pipe or sphere sizes for input to the network model: Values of the coefficients a and b for pipes and spheres are given in Table 1.The linear cumulative probability distribution of Eq.6 corresponds to a uniform probability distribution between a minimum radius of b/a and a maximum radius of (1 + b)/a , and a mean radius of (1 + 2b)/(2a).With the values of a and b presented in Table 1, the pipe radii varied between 1.35 e-9 m and 5.00 e-5 m, the sphere radii varied between 1.35 e-9 m and 3.44 e-4 m and the mean sphere radius r sm was 1.72 e-4 m.
Having defined the mean sphere radius r sm , the values of l SVE (the SVE side length) and d min (the minimum distance between the randomly positioned points forming the Delaunay tessellation) were selected to give an appropriate number m of spheres within the SVE and to match the reported porosity of the sandstone of 0.21.The values of l SVE and d min shown in Table 1 produced a SVE containing 1560 spheres.A series of 100 analyses was performed.The same network geometry was used for all 100 analyses, but different sphere and pipe radii generation and allocation were undertaken for each analysis.The number of the seeding spheres used for the initiation procedure described in Section 3.2 was 1% of the total number of spheres, so that drying and wetting curves were simulated over a range of S r from 0.99 to 0.01.The same increments of suction were applied for all 100 analyses.

Results
The circular data points in figure 5 show the results of the network model simulations in S r -P c space.Values of S r represent the average from the 100 analyses.6 show the results of the network model simulations re-plotted in the space of S r against back-calculated radius r (where Eq. 4 was used to convert values of P c to corresponding values of r). Figure 6 provides insights into some of the mismatches between model simulations and experimental results in figure 5.
In figure 6, comparison of the network model simulations of S r against back-calculated r with the cumulative probability distributions of spheres and pipes used as input for the network model (shown by the straight lines in figure 6) highlights the influence of network connectivity on wetting and drying curves.For example, the network model simulation of the wetting curve predicts lower values of S r than those corresponding to the cumulative probability distribution of sphere sizes, because network connectivity means that some smaller spheres cannot be intruded with the wetting fluid at their expected value of P c , because access of the wetting fluid to these spheres is shielded by some larger spheres.As the degree of saturation increases, the network model simulation converges with the cumulative probability distribution of spheres, because this shielding effect reduces as more spheres are filled with the wetting fluid.
The model simulations in figure 5 could be tuned to better match the experimental results by adjusting the values of some of the input parameters in Table 1.This would mean, however, that the method had reverted to no more than a curve fitting exercise.Proper assessment of the achievements and limitations of network model simulations of retention behaviour can be achieved only by using inputs for the sphere and pipe probability distribution functions derived from independent imaging of the material microstructure, rather than from backanalysis of retention behaviour from mercury porosimetry tests.
One significant point to emerge from figure 5 is that the match achieved between model simulations and experimental results was better than in previous network model simulations [10] which did not use the pore volume scaling technique, even though these earlier simulations included substantial tuning of input parameter values to try to achieve the best possible match.This provides encouraging evidence on the usefulness of the pore volume scaling technique.

Conclusions
A network model was used to study retention behaviour in unsaturated materials.Original aspects of the modelling included a form of periodic boundary conditions applied for the first time to a 3-D transport network, a new pore volume scaling technique and a new procedure for initiating drying or wetting paths in a network model employing periodic boundary conditions.The model is suitable for future extension to hydromechanical and multi-scale analysis.The model was applied to the simulation of a mercury porosimetry test on sandstone, with linear cumulative probability distribution functions selected for both sphere and pipe sizes, based on back-analysis of the experimental test results.Comparison of model simulations with both experimental results and the cumulative probability distributions used for sphere and pipe sizes provided useful insights into the roles of network connectivity on drying and wetting retention curves predicted by network models.

Figure 1 .
Figure 1.Schematic presentation of four Voronoi polyhedra and the corresponding Delaunay tetrahedron.The spheres and pipes corresponding to a Voronoi polyhedron face are also presented.

Figure 2 .
Figure 2. Example of a 3-D fully periodic network where the grey and the black edges represent the pipes and the SVE edges, respectively.
of the SVE and therefore the porosity n of a cubical SVE of side length l SVE containing a total of m spheres is simply given by

Figure 3 .
Figure 3.A simple 2-D network subjected to drying.The hatched areas represent the wetting fluid.The drying fluid is allowed to enter the network from pipe 1.

Figure 4 .
Figure 4.A simple 2-D network subjected to wetting.The hatched areas represent the wetting fluid.The wetting fluid is allowed to enter the network from pipe 1.

Figure 5 .
Figure 5.Comparison of the experimental results ([5]) to the mean numerical results in S r -P c space.

Figure 6 .
Figure 6.Comparison of the experimental results ([5]) and the cumulative radii distribution curves to the mean numerical results in the S r -r space.
Comparison with the experimental results shows that the network model simulation of the drying curve overestimates values of S r at low values of P c and underestimates values of S r at high values of P c .The network model simulation of the wetting curve underestimates values of S r at all but the lowest experimental value of P c .The circular data points in figure