Capillary pressure–saturation relationships for porous granular materials: Pore morphology method vs. pore unit assembly method

In studies of two-phase ﬂow in complex porous media it is often desirable to have an estimation of the capillary pressure–saturation curve prior to measurements. Therefore, we compare in this research the capability of three pore-scale approaches in reproducing experimentally measured capillary pressure– saturation curves. To do so, we have generated 12 packings of spheres that are representative of four different glass-bead packings and eight different sand packings, for which we have found experimental data on the capillary pressure–saturation curve in the literature. In generating the packings, we matched the particle size distributions and porosity values of the granular materials. We have used three different pore-scale approaches for generating the capillary pressure–saturation curves of each packing: i) the Pore Unit Assembly (PUA) method in combination with the Mayer and Stowe–Princen (MS–P) approximation for estimating the entry pressures of pore throats, ii) the PUA method in combination with the hemisphere approximation, and iii) the Pore Morphology Method (PMM) in combination with the hemisphere approximation. The three approaches were also used to produce capillary pressure–saturation curves for the coating layer of paper, used in inkjet printing. Curves for such layers are extremely diﬃcult to deter- mine experimentally, due to their very small thickness and the presence of extremely small pores (less than one micrometer in size). Results indicate that the PMM and PUA-hemisphere method give similar capillary pressure–saturation curves, because both methods rely on a hemisphere to represent the air– water interface. The ability of the hemisphere approximation and the MS–P approximation to reproduce correct capillary pressure seems to depend on the type of particle size distribution, with the hemisphere approximation working well for narrowly distributed granular materials. © 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The relationship between capillary pressure and saturation plays an important role in many applications of porous materials, such as inkjet printing ( Rosenholm, 2015 ), water flow in hygienic products ( Diersch et al., 2010 ), water distribution in fuel cells ( Qin, 2015 ), oil and gas reservoir engineering, and unsaturated soils ( Lins and Schanz, 2005 ). This relationship depends on both the geometry and dimensions of the pore space ( Likos and Jaafar, 2013;Mousavi and Bryant, 2012;Øren and Bakke, 2003;Torskaya et al., 2014 ).
Usually, it is possible to represent the pore space of a granular material by a network of pore bodies and pore throats. The pore bodies account for most of the void volume, while the pore throats, which connect neighbouring pore bodies, act as narrow pathways that account for the resistance to flow. The network of imaging. Using image analysis techniques, a 3-dimensional highresolution reconstruction of the pore space can be obtained in binary mode. The resolution of this method is in the order of nanometer. However, the acquisition of FIB-SEM images is time consuming (between 4 and 12 h depending on sample size and slicing thickness) and the analyzed volume is typically very small ( Aslannejad et al., 2017;Bultreys et al., 2016 ). An example of a non-destructive method is X-ray computed tomography that provides 3-dimensional structural information, with a resolution in the range of micrometer to millimeter scale. It can visualize not only the pore space but also fluid distributions inside the porous material (e.g. Culligan et al., 2004 ). However, direct measurements of the pore geometry are not trivial and it remains time consuming, especially for visualization of numerous samples. As an alternative, for the purpose of studying various factors that affect hydraulic properties of granular materials, one may use numericallygenerated packings of spheres.
One way to generate a packing is to place spheres randomly in a stable configuration by trial and error until a desired porosity is obtained. Bakke and Øren (1997) generated packings of randomly located spheres, before compaction was applied by scaling all particle locations in the vertical direction. The pore geometry was extracted using image analysis techniques on the packing of spheres. Based on the resulting pore geometry, a pore network was constructed that was used to obtain capillary pressure-saturation curves and the relative permeability curves. Mousavi and Bryant (2012) generated a random packing of ductile spheres in which cementation was modelled by increasing the radii of the spheres while accounting for a net transport of solidphase, due to pressure-dissolution, from finer grains to coarser grains. Once a packing was produced, triangulation was applied to extract the pore geometry and subsequently to construct the capillary pressure-saturation curve. They studied the effect of cementation on the capillary pressure-saturation relationship.
Alternatively, one may start with a collection of sparsely placed spheres and use the Discrete Element Method (DEM) to determine their stable configuration. DEM is based on determining the force interactions between different spherical particles during their displacement and compaction. DEM has been used to simulate geo mechanical processes, such as triaxial tests ( Belheine et al., 20 09;Plassiard et al., 20 09;Widuli ński et al., 2009 ), soil tillage ( Shmulevich, 2010 ), as well as industrial processes such as particle flows in silos ( Langston et al., 1995 ). DEM has been used to construct particle packings with a predefined particle size distribution and porosity value ( Tong et al., 2012;Yuan et al., 2016 ). One disadvantage of DEM is that the shape of individual particles is simplified to a sphere, which may not be appropriate for irregularly shaped particles. Nevertheless, packings of spheres have proved to be valuable for the study of processes in granular materials.
The pore space in spherical particle packings is often extracted using Delaunay tessellation, also known as triangulation, where centres of four neighbouring particles are connected to form a tetrahedron. This tetrahedron is then assumed to enclose one pore unit such that an assembly of pore units represents the whole pore space ( Mason and Mellor, 1995 ). We refer to this method as the Pore Unit Assembly (PUA) method. Gladkikh and Bryant (2005) have used triangulation to extract the pore space of a Finney packing, which is a packing of equallysized spheres whose locations have been measured ( Finney, 1970 ), to obtain the capillary pressure-saturation curve for imbibition. Yuan et al. (2016) generated packings of spheres, using DEM, which represented the glass-bead packing used in experiments of Culligan et al. (2004) . Then, they used PUA to construct the drainage capillary pressure-saturation curve for that packing. They found good agreement with experimental data of Culligan et al. (2004) . The concept of PUA can be gener-alized such that it can be applied to other porous materials types.
Another method to compute the capillary pressure-saturation curve from both direct pore geometry measurements and packings of spheres is the Pore Morphology Method (PMM). In PMM, the air-water interface is represented by a hemisphere where the radius of that sphere is related to the capillary pressure, using Young's Laplace equation. PMM divides the pore structure into voxels. Then, for a certain capillary pressure, PMM identifies where an air-water interface would fit in the pore space while maintaining the connectivity of the air phase to its reservoir. At each capillary pressure, the distribution of water and air phase is known, such that the capillary pressure-saturation curve can be calculated ( Hilpert and Miller, 2001 ). Vogel et al. (2005) have compared the capillary pressure-saturation relationship of PMM with these obtained from a pore-network model and Lattice Boltzmann simulations, where they solved for two-phase Navier-Stokes on the actual pore geometry. The pore network contained a 3-dimensional network of capillary tubes with the same statistical information on pore throat size distribution and coordination number as used in PMM computations. In that research, a 3-dimensional reconstruction was made of a sintered glass sample using X-ray tomography results. The pore space was then used to compute the capillary pressure-saturation curve using the three different models. Results showed that Lattice Boltzmann gives lower values of capillary pressures compared to that of the PMM and the pore network model. The capillary pressure-saturation curves for PMM and the pore network model were approximately the same, because the pore network model inherited the pore structure that was extracted in the pore morphology method. However, the question remains on how well both models can predict the capillary pressuresaturation curve of granular materials if the pore network is not directly based on the statistical information from PMM, but rather on the discretization of the pore space into pore throats and bodies.
An essential factor in determining the capillary pressuresaturation curve is the shape and curvature of the air-water interface ( Vidales et al., 1998 ). Air-water interfaces in granular materials can have complex shapes, whose curvatures are difficult to determine. In pore-scale models, the curvature, at which air can invade a water-saturated pore throat, can be estimated by either the Mayer and Stowe-Princen (MS-P) approximation or by the so-called hemisphere approximation. The MS-P approximation assumes that the main terminal meniscus (MTM) has just passed through a pore throat opening (so, it is inside the pore throat) and corners in the invaded part of the pore throat are still filled by water. Then, the capillary pressure of the MTM can be assumed equal to that of the arc menisci (AM) which are present in the corners. The radius of the curvature of AM in the cross-section of a pore throat can be determined by a force balance. The second curvature, along the pore throat, is assumed to be flat. This assumption is valid only if the pore throat is elongated. Nevertheless, Prodanovi ć and Bryant (2006) have used the level-set method to study the shape of an air-water interface in pore throats formed between spherical particles. They showed that the 3-dimensional curvature of air-water interface in a throat can be obtained by the MS-P approximation. The hemisphere approximation assumes that the airwater interface has a spherical shape, and thus assumes a contact angle of zero. In PMM, the air-water interface is represented by a hemisphere ( Hilpert and Miller, 2001 ) while in PUA method and pore-network models both the hemisphere and the MS-P approximations have been used.
To investigate the ability of PUA and PMM to estimate capillary pressure-saturation curves of granular materials, we studied 12 granular materials for which experimental data are reported in the literature. In addition, we studied the coating layer of paper which has a mean pore size of 180 nm ( Aslannejad et al., 2017 ). We generated for each granular material a spherical particle packing by matching the particle size distribution and the porosity value of the granular material of interest, using DEM. This packing was then imported into both PMM and PUA to construct capillary pressure-saturation curves for drainage using three different methods: (i) using PUA where entry pressures are given by the Mayer and Stowe-Princen method, (ii) using PUA where entry pressures are given by the hemisphere approximation, and (iii) using PMM in combination with the hemisphere approximation.
The aims of this research were as follows: a) determine how capillary pressure-saturation curves from PUA and PMM compare to those from experiments and how this is affected by the choice of approximation for the entry pressure, using either the MS-P approximation or hemisphere approximation; b) Determine how the shape of the air-water interface is affected by the particle size distribution. To do this, the following three comparisons were made: (i) we compared capillary pressure-saturation curves constructed by PMM to those of PUA, both of them using the hemisphere approximation; (ii) we compared the hemisphere approximation with the MS-P approximation, both used in PUA, to understand what their main differences are, and (iii) we compared capillary pressure-saturation curves obtained with PUA, using both MS-P and the hemisphere approximations, to existing experimental data to determine their ability of reproducing experimental results.

Discrete element method
To construct packings of spheres with a predefined porosity and particle size distributions, we use open source software Yade-DEM ( Šmilauer et al., 2015 ), which is a 3-dimensional particle model based on the Discrete Element Method (DEM). In DEM, individual particles are represented by spherical elements. For each discrete element, a force and rotation balance is computed which is based on interactions with boundaries and other discrete elements. At the contact between two discrete elements, normal displacement, tangential displacement, and sliding can occur. To explain this, consider two particles that are in contact and pushed towards each other due to interactions with their surrounding particles. Due to the compression at their contact, both particles will flatten at their contact. As a consequence, an elastic force arises because of the tendency of both particles to restore their initial shape. This elastic force is simulated as a spring while the discrete elements act as dashpots and therefore the model is often referred to as a springdashpot model. Furthermore, tangential displacement may occur at the contact, which can give rise to sliding, depending on the static friction angle of the material. For more information on DEM, the reader is referred to Šmilauer et al. (2015) .

Pore unit assembly method 2.2.1. Pore unit geometry
The pore unit assembly method (PUA) that we use in this research is a module in Yade-DEM Sweijen et al., 2016;Yuan et al., 2016 ). PUA relies on regular triangulation to subdivide the pore space of packings of spheres into pore units, see Fig 1 . In a triangulation, the centres of four neighbouring particles form the vertices of a tetrahedron, see Fig 1 b. The void space that is present within a tetrahedron is referred to as a pore unit . Each tetrahedron is connected to four other tetrahedra. The facet that is shared between two pore units of adjacent tetrahedra represents the narrowest opening between two tetrahedra, which we refer to as a pore throat . A pore throat has zero length but has a cross-sectional area.

Construction of capillary pressure-saturation curves
To construct capillary pressure-saturation curves for drainage, it is essential to know the entry pressure that is required for air to invade a water-saturated pore unit via one of its four pore throats. Consider an air-saturated pore unit i at an air pressure P i and a water saturated pore unit j at a water pressure P j . For air to invade pore unit j , the pressure difference between pore units i and j ( P i − P j ) should be larger than the entry pressure of pore throat ij , namely P e i j . To compute the entry pressure of pore throat ij , two different methods can be used. Both methods are based on the Young-Laplace equation, which relates the entry pressure to the curvature of the air-water interface: where r and r are two principal radii of curvature of the airwater interface and γ aw is the air-water interfacial tension.
The first method is known as the Mayer and Stowe-Princen method ( Mayer and Stowe, 1965;Princen, 1969 ). We refer to this method as PUA-MS-P . To explain this method, consider an elongated pore throat ij , between air-saturated pore unit i and water saturated pore unit j (see Fig 2 a). Let's assume that air has just invaded it. The entry pressure can be computed from the balance of forces acting on the air-water interface in the cross-section of pore throat ij , by assuming that the capillary pressure just after invasion is the same as the entry pressure of that pore throat (see e.g. Ma et al., 1996 ). The radius of curvature of the arc menisci (in the pore throat corners) R MSP i j is obtained based on the pore throat that was found during triangulation. The particle locations of the three spheres that form a pore throat determine the shape of that pore throat and thus determine the value of R MSP i j (see Appendix A ).
Note that radius R MSP i j is supposed to be smaller than the radius of the inscribed circle of a pore throat R in i j ( Yuan et al., 2016 ). This, however, may not occur when particles do not touch, as illustrated in Fig 2 b. In the MS-P approach it is assumed that the pore throat is elongated such that one principal radius in Eq. (1) goes to infinity. Thus, Eq. (1) can be approximated by: For more details on the implementation of this method into DEM, the reader is referred to Yuan et al. (2016) .
The second method to compute the entry pressure of a pore throat is to assume that the air-water interface has the shape of a hemisphere, with a contact angle of zero, see Fig 2 c. We refer to this method as PUA-hemisphere . The radius of the hemisphere corresponds to the inscribed circle in a pore throat R in i j . Then, the principal radii in Eq. (1) are assumed to be both equal to R in i j such that Eq. (1) can be recast into: For both MS-P approximation and hemisphere approximation, the model by Sweijen et al. (2016) was employed to obtain the capillary pressure-saturation curve for drainage. In that model, the saturation was assumed to be binary; i.e. a pore unit is assumed to be saturated either with water or air. Thus, pendular rings and water films are neglected. The simulation procedure to obtain capillary pressure-saturation curves was as follows. Initially, all pore units were considered to be saturated with water, except for the pore units along the air reservoir which were saturated with air. Next, a capillary pressure was imposed on the packing. This was done by setting the pressure difference between the air and water reservoirs ( P air − P w ) equal to the value of that capillary pressure. Then, for each pore throat in between an air saturated pore unit Illustration of a pore unit: a) four spheres enclose one pore unit, b) a tetrahedron that is a result of regular triangulation, with pore throat radius R ij and c) illustration of the inscribed sphere of a pore unit with a radius R i .

Fig 2.
Schematic representation of air-water menisci in two different approximations for the estimation of the entry pressure of a pore throat; a) the Mayer and Stowe-Princen approximation, where the air-water interface has a complex shape in an elongated pore throat, b) an example of a pore throat with non-touching particles, and b) the hemisphere approach where the air-water interface is assumed to be spherical.
i and a water saturated pore unit j , the following criterion was evaluated: if the imposed capillary pressure was larger than the entry pressure ( P e i j ), air would invade pore unit j . If invasion occurred, then the criterion was checked again for all pore throats connecting to pore unit j . This was continued until no air invasion would occur. The saturation corresponding to that imposed capillary pressure was measured. Then, a new capillary pressure was imposed, the model checked again for invasions, and the new saturation was measured. This resulted in a capillary pressuresaturation curve. The algorithm does not allow for displacement of disconnected water-filled pore units, resulting in a residual water saturation after drainage. For more details, the reader is referred to Sweijen et al. (2016) .

Pore morphology method
The third method that we use to calculate capillary pressuresaturation curves is the Pore Morphology Method to which we refer to as PMM . For this we used GeoDict software. The pore morphology method determines the stationary air and water distribution in a porous material for a given capillary pressure. In static conditions, the capillary pressure is the same as the pressure difference between the air reservoir and the water reservoir, which are at a pressure P air and P w , respectively. To determine the pore space that is accessible to air (during drainage) the Young-Laplace equation is employed, where the air-water interfacial area is as-sumed to be a hemisphere with radius R c : To identify the locations in the pore space where a hemisphere with radius R c would fit, an erosion algorithm is employed, which is normally used in image analysis procedures. Considering the connectivity of air to the air-reservoir, the capillary pressure-saturation curves were constructed. For more information, the reader is referred to Hilpert and Miller (2001) .

Simulation methods
In this research, we study granular materials that are reported in the literature, for which capillary pressure-saturation curves have been measured. For each granular material, we generated a packing of spheres that matched its porosity value and its particle size distribution. For each packing, capillary pressure-saturation curves were produced using the three different methods as described previously.

Materials
The particle size distributions of the 12 materials are listed in Table 1 . We have chosen to simulate 4 glass bead samples and 8 sand samples. Out of 8 sands, four were Ottawa sands (samples OS12/20, OS20/30, OS30/40 and OS40/50). They had different particle size distributions, at a constant porosity value of 0.34  ( Schroth et al., 1996 ), which enabled us to only study the effect of the particle size distribution. In addition, two sands that are reported in Oostrom et al. (1999) were studied, as they also had the same porosity value of 0.34 and originated from the same sand type but had different filter sizes (12S and 14S). The two remaining sands (F75, RS) had higher porosity values and different particle size distributions. They were used to study the applicability of the models on different types of sands. Two samples of glassbeads were chosen because they have been used in previous studies (GB1 and GB2), which enabled us to cross-verify our results ( Culligan et al., 2004;Hilpert and Miller, 2001 ). The other two glass bead samples (GB3 and GB4) were chosen because they had varying particle sizes and varying porosity values. Most granular materials had a porosity value between 0.34 and 0.36, but 4 granular materials had a porosity value higher than 0.38. All experimental capillary pressure-saturation curves were scaled to air-water systems using a surface tension of 0.072 Nm −1 and a contact angle of zero degrees.
Finally, we generated a packing of spheres that corresponds to the coating layer of paper, which is reported in Aslannejad et al. (2017) . They extracted the pore geometry using FIB-SEM imaging and they subsequently constructed the capillary pressuresaturation curve using PMM directly on the extracted pore geometry. We use their reported values of Van Genuchten parameters to evaluate the approach.

Generation of a pack of spheres
In our DEM simulations, we have used the linear contact mechanics, as introduced by Cundall and Strack (1979) . The mechanical parameters were adopted from Belheine et al. (2009) , who have fitted DEM simulation of various triaxial tests to experimental results. The particle stiffness was set to 9.8 × 10 8 N m −1 , the ratio between the normal and tangential force to 0.04, and the density was set to 2600 kg m −3 .
First, a cloud of 40 0 0 spheres was generated in a modelling domain (0.02 × 0.02 × 0.02 m ³) that was larger than the final packing dimensions, which was typically in the order of 1 mm 3 to 1 cm 3 . The particle sizes were chosen from the particle size distribution of the granular material of interest. Then, an arbitrarily large friction angle (50 °) was assigned to the particles. Next, a confining pressure was applied to all 6 boundaries of the box such that the cloud of particles got compacted to a stable configuration, which had a porosity value of approximately 0.43. Then, the friction angle was decreased in small steps such that the particles were gradually rearranged into a tighter packing until the desired porosity was achieved.

Capillary pressure-saturation relationship
The packing of spheres was imported in PMM whereas PUA was directly employed within Yade-DEM. During capillary pressuresaturation calculations in all methods, one boundary of the modelling domain of DEM was assumed to be an air reservoir at an air pressure P air , and the opposite boundary was assumed to be a water reservoir at a water pressure of P w . The air pressure was increased such that air could invade increasingly smaller pores.

Van Genuchten function
To allow for quantitative comparison of all capillary pressuresaturation data, we have fitted each curve with the Van Genuchten function ( Van Genuchten, 1980 ): in which α [Pa −1 ] is proportional to the inverse of the entry pressure, n is related to the pore size distribution, and S e is the effective saturation, which is defined for primary drainage as: Here, S is the saturation, S r is the irreducible water saturation after drainage.

Verifying the packing procedure
To verify our approach of constructing capillary pressuresaturation curves of packings of spheres, we have reproduced a capillary pressure-saturation curve for the packing of spheres that corresponds to the glass beads used in Hilpert and Miller (2001) , namely sample GB1. They measured the primary drainage curve and simulated it using PMM. Both their experimental and simulated curves (for 800 ³ voxels), as well as our PMM simulated curve are shown in Fig 3 . The capillary pressure-saturation curves that were constructed using PMM in Hilpert and Miller (2001) and in our simulations are similar, indicating that the packings of spheres are similar. Albeit we used DEM to make a packing of spheres while they used a random insertion algorithm that was developed by Yang et al. (1996) . However, none of the simulated curves can match the experimentally measured curve. Both PMM curves have a smooth transition at the entry pressure, while experiments show a discontinuity of the curve at the entry pressure. In other words, the simulated curves can be fitted best Capillary pressure-saturation data for a glass-bead packing (GB1) together with curves constructed using the pore morphology method (PMM). Simulated curves by Hilpert and Miller (2001) were for 800 3 voxels. with the Van Genuchten formula, because of the continuous behaviour of the curve around the entry pressure, whereas the experimental behaviour can be fitted best by the Brooks and Corey formula ( Brooks and Corey, 1964 ), because of its discontinuity around the entry pressure. In the saturation range of 0.2 to 0.8 the capillary pressure of the simulated curves corresponds reasonably well to the experimental data. However, at low saturations ( < 0.20), the capillary pressure is overestimated by the PMM, because the hemisphere approximation does not hold for pendular rings ( Hilpert and Miller, 2001 ). For pendular rings, the two principal radii r and r in Eq. (1) are not equal; in fact, one of them is even negative ( Scholtès et al., 2011 ). Fig 4 shows the capillary pressure-saturation curves by PMM, PUA-MS-P and PUA-hemisphere, for the same packing as in Fig 3 . PUA-hemisphere results in a curve that fits the experimental data better than the curve by PMM, especially for saturations larger than 0.8. Nonetheless, the PMM curve is still relatively close to that of PUA-hemisphere, because both methods assume the air-water interface to have a spherical shape.

Hemisphere approximation: Pore morphology method vs. pore unit assembly model
As we discussed in Section 4.1 , the capillary pressure-saturation curves are close to each other for PMM and PUA-hemisphere simulations on glass bead packing GB1. To confirm this for all packings, we have plotted the values of Van Genuchten parameter α for all 12 granular materials obtained with PMM against values obtained with PUA-hemisphere (see Fig 5 a). In addition, the Van Genuchten parameters are reported in Table 2 . Indeed, both mod-els give approximately the same value of α and thus give approximately the same value of entry pressure. However, the values of Van Genuchten parameter n obtained from the two models are not the same, as can be seen from Fig 5 b. The PUA-hemisphere method gives flatter capillary pressure-saturation curves, and thus larger values of n , than the PMM.

Hemisphere approximation vs. Mayer and Stowe-Princen approximation
To study the capillary pressure-saturation curves obtained from PUA-hemisphere and PUA-MS-P approximations, we show in Fig 6  both curves for 8 granular materials having porosity values between 0.34 and 0.36. The results show that the hemisphere approximation gives an entry pressure that is twice as large as that of the MS-P approximation. To identify the underlying cause of this difference, we study the distributions in R ins i j and R MS−P i j for sand OS12/20, see Fig 7 . Typically, the radii of curvature for arc menisci, R MS−P i j , are smaller than that of R ins i j ( Bico and Quéré, 2002 ), which is also the case for throats with a triangular crosssection ( Mason and Morrow, 1991 ). This is reflected in Fig 7 , where MS-P results in an abundancy of small values of R MS−P i j ( ±44 μm) that is smaller than the smallest value of R ins i j , where R ins i j has an abundancy at 79 μm. However, both methods have a similar distribution for larger pore throats ( > 200 μm), which is related to the triangulation of packings of spheres. One throat is enclosed by three spheres that either touch each other or not, depending on the force consideration in DEM. Therefore, some pore throats are not perfectly enclosed by solid and thus finding the value of R MS−P i j is not possible as the distance between two spheres may be too large to have R MS−P i j smaller than R ins i j . This occurs in 8.5% of the pore throats, which are typically relatively larger pore throat ( > 150 μm) Therefore, when the difference between air pressure and water pressure just exceeds the entry pressure of a sample, air will invade larger pore throats that have values of R MS−P i j , which are close to that of R ins i j . Therefore, the only difference between the MS-P and hemisphere approximation is the evaluation of entry pressure (using either Eqs. (2) or (3) , respectively), which indeed is a factor of two; this of course is an artefact of simulations rather than a physical phenomenon. Thus, we may conclude that MS-P cannot always be applied for simulations in the Pore Unit Assembly.

Comparison of simulated curves and experimental data
Fig 6 also shows that the drainage curves of PUA simulations, both the hemisphere approximation and the MS-P approximation, lack the steepening of the curve near residual water saturation, which are present in experimentally measured curves. This is mainly related to model simplifications; the particles are spherical and flow of water along corners is not simulated. Indeed, the pore throat size distribution of packings used in experiments may differ from that in simulations because particles will have a different shape than spheres. Also, corner flow is responsible for additional drainage which is not considered in the simulations and, therefore, the residual water saturation is slightly overestimated by the model.
Results shown in Fig 6 indicate that PUA-hemisphere method approximates well the entry pressure of five granular materials (OS12/20; OS20/30; OS30/40; OS40/50; GB1), but gives poor approximations for 12S, 14S and GB2. However, the MS-P approximation reproduces the capillary pressure-saturation curve for glassbead sample GB2 satisfactory. Glass bead sample GB2 has been previously investigated by Yuan et al. (2016) , who have constructed  Van Genuchten parameters for all 12 granular materials and the coating layer (Van Genuchten parameters for experiments have been extracted from the literature, they are compared to parameter values fitted to capillary pressure-saturation curves obtained from simulations using three different methods).

Type
Name Experiments PMM-Hemisphere PUA -Hemisphere PUA -MS-P  Brooks and Corey (1964) , S e = ( P e P c ) λ , gave a better fit with best-fit values being: P e = 4312 Pa and λ= 10.7. * * In the original work by Oostrom et al. (1999) the Brooks and Corey (1964) parameters were reported for tricholoethylene and water, which was converted to air and water by correcting for interfacial tensions (0.035 Nm −1 to 0.072 Nm −1 , respectively). The scaled-version of the Brooks and Corey parameters were then converted to Van Genuchten parameters. * * * Van Genuchten parameters values were obtained using Pore Morphology Method and based on directly measured pore geometry and we corrected the reported values for a contact angle of 45 °i ts capillary pressure-saturation curve using PUA-MS-P simulations. Indeed, they also found a good match between experiments and simulations, similar to our results. Moreover, Joekar-Niasar et al. (2010) have also successfully reconstructed the capillary pressure-saturation curve of GB2 using a pore-network model, based on a pore network that was obtained from X-ray tomography. They also used the MS-P approximation.  9 ). One would expect a linear correlation between α and D50 ( Benson et al., 2014 ). Indeed, Fig 9 shows that all data points fall on a straight line except for the three samples. We be-lieve that the deviation of these samples from the line is not an experimental error.
The fact that the three samples behave differently in experiments implies that the nature of the granular material could affect the relation between α and D50; of course, the type of particle size distribution is another factor . Fig 10 a shows the particle size distributions of all Ottawa sands for which the hemisphere approximation works well. These Ottawa sands have a relatively narrow and normal distribution of particle sizes. However , Fig 10 b shows the wider distributions of samples 12S and 14S for which the MS-P approximation yields better results than the hemisphere approximation. Finally, Fig 10 c shows the very irregular particle size distributions of sample GB2 for which the MS-P approximation works well. We may postulate that the particle size distribution may affect the shape of pore throats and thus the shape of air-water interfaces. For example, elongated pore throats may exist in packings that have a wider particle size distributions. In elongated pore throats, the shape of the air-water interfaces may substantially deviate from that of a hemisphere, such that the MS-P approximation may be more applicable. But, for narrow distributions where the Capillary pressure-saturation curves for the granular materials that have a porosity value around 0.35. The dots represent experimental data, the lines represent simulations with PUA-hemisphere method (solid line) and PUA-MS-P method (dashed line). Note that for sands 12S and 14S, dots represent the fitting of Van Genuchten formula from data reported in Table 2 and Oostrom et al. (1999) .  pore throats are thin, the hemisphere approximation seems work well for simulations in PUA. However, another explanation for the different behaviour of the three samples in our modelling effort s is that the packing of spheres does not represent the real granular material well, leading to a wrongly estimated capillary pressure-saturation curve. But, our results on glass-beads GB1 and GB2, which are best reproduced by the hemisphere approximation and the MS-P approximation, respectively, have been cross-referenced with other simulations and experiments, implying that our simulations are in accordance with other studies. Another explanation is that we have chosen samples that, in hindsight, tend to be better represented by the hemisphere approximation. While, for another set of samples the hemisphere approximation may not work well; it remains a process of sampling from a large population of granular materials. Therefore, in the scope of this work, it is not possible to identify with certainty why the three samples behave differently from the others.
To exemplify the concept of generating a capillary pressuresaturation curve prior to experiments, we evaluate the coating layer of paper, which is a typical example of cases where direct measurements of the capillary pressure-saturation curves are not feasible. Moreover, grains in the coating layer are irregularly shaped ( Aslannejad et al., 2017 ) and thus the representation of  Table 2 (last entry) was actually obtained by Aslannejad et al. (2017) from PMM simulations on the real pore geometry of the coating layer, obtained from FIB-SEM measurements. We also determined the value of α, for a packing of spheres with the same particle size distribution, using PMM, PUA-MS-P, PUA-hemisphere methods. We found that the α value was underestimated by the hemisphere approximation while overestimated by PUA-MSP method. The fact that PMM was not capable to reproduce the α value from Aslannejad et al. (2017) is caused by the simple pore geometry used in our results; i.e. that of a packing of spheres. Nevertheless, the value of α from the direct pore geometry is in between the α value by PUA-MS-P and PUA-hemisphere. Of course, the prediction of α values is better for granular materials with more spherical particles (i.e. glass-beads and filter sands) than for irregularly shaped particles.

Conclusion
In this research, we have constructed capillary pressuresaturation curves for drainage of spherical particle packings using three different models with varying assumptions on the shape of the air-water interface, namely: (i) the Pore Unit Assembly method (PUA) in combination with the Mayer and Stowe-Princen method, (ii) PUA in combination with the hemisphere approximation or (iii) the Pore Morphology Method (PMM) in combination with the hemisphere approximation. We have generated spherical particle packings for 12 granular materials from the literature and the coating layer of a paper. From each packing, the capillary pressuresaturation curve was constructed using the three methods and the results were compared to experimental data from the literature. Results show that PMM and PUA-hemisphere give similar capillary pressure-saturation curves, albeit that PUA gives flatter curves than PMM. In these simulations, the hemisphere approximation results in an entry pressure value twice that of MS-P. Comparison with experimental data from the literature reveals that for five out of eight samples the hemisphere approximation reproduces capillary pressures well. The five samples have a narrow particle size distribution and show a different relation of experimentally determined α with D50 than the other three samples, which have a relative wide particle size distribution. Therefore, we suggest that the type of particle size distribution affects the type of invasion criterion that should be employed in capillary pressure-saturation computations, in the pore unit assembly method.

Acknowledgment
The first author gratefully acknowledges financial support from the Technology Foundation STW, the technological branch of the Netherlands Organisation of Scientific Research, NWO, and the Dutch ministry of Economic Affairs under contract no. 12538, entitled Interfacial effects in ionized media.
The second and third authors would like to thank European Research Council (ERC) for the support they have received under the ERC Grant Agreement no. 341225 . Constructive comments by three anonymous reviewers helped to improve the manuscript significantly.

Appendix A -Mayer and Stowe -Princen approximation
The radius of curvature of arc menisci in the cross-section of a pore throat is computed using the MS-P approximation. Fig 2 a shows a cross-section where a main terminal meniscus has just entered the pore throat and corners in the invaded part of the pore throat are still filled by water, forming arc menisci. The following quantities can be identified in Fig 2 a: the length of air-water contact line ( L aw ), the length of the air solid contact line ( L as ), the cross-sectional area that is filled with air ( A n i j ) and the radii of the air-water interface in the plane of the pore throat ( R MSP i j ). The values of A n i j , L aw and L as are found on the pore geometry resulting from triangulation of the packing of spheres, for the computation of these values, see Yuan et al. (2016) . For computing R MSP i j , we note that surface tensions should in balance with the pressure difference between pore unit i and j , such that: P i − P j A n i j = L aw γ aw + L as γ as − L as γ ws