{'fraction_non_alphanumeric': 0.029201900851363188, 'fraction_numerical': 0.012820696883120613, 'mean_word_length': 4.660165772793759, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ability of a pore network model to predict fluid flow and drag in saturated granular materials

The local flow field and seepage induced drag obtained from Pore Network Models (PNM) is compared to Immersed Boundary Method (IBM) simulations, for a range of linear graded and bimodal samples. PNM were generated using a weighted Delaunay Tessellation (DT), along with the Modified Delaunay Tessellation (MDT) which considers the merging of tetrahedral Delaunay cells. The local pressure field was very accurately captured in all linear graded and bimodal samples. Local flux (flow rate) exhibited more scatter, but the PNM based on the MDT clearly provided a better correlation with the IBM. There was close similarity in the network shortest paths obtained from PNM and IBM, indicating that the PNM captures dominant flow channels. Further, by overlaying the PNM on a streamline profile, it was demonstrated that local pressure drops coincided with the pore constrictions. A rigorous validation was undertaken for the drag force calculated from the PNM by comparing with analytical solutions for ordered array of spheres. This method was subsequently applied to all linear graded and bimodal samples, and the calculated force was compared with the IBM data. Linear graded samples were able to calculate the force with reasonable accuracy, while the bimodal sample exhibited slightly more scatter.


Introduction
A detailed understanding of the fluid flow at the level of the individual particles and the drag induced by the seepage flow in saturated sands is important in applications including sand production in hydrocarbon extraction (Climent et al., 2014), filter design (Huang et al., 2014), piping erosion (Tao and Tao, 2017) and suffusion (Kawano et al., 2018;Hosn et al., 2018) in dams and flood embankments and fundamental studies of liquefaction (El Shamy et al., 2010). The topology of the void space between the particles can now be determined relatively accurately using micro-computed tomography. Moreover, model soils can be developed using Discrete Element Method (DEM) simulations, from which particle position data and contact locations are readily available. This topology data can be used to develop fully resolved flow simulations using Computational Fluid Dynamics (CFD), in which the fluid motion is described at a scale that is smaller than the individual particles, with particle surfaces considered as boundaries to the fluid flow. In these simulations, the Navier-Stokes equations can be solved using the finite element method (Glowinski et al., 2001), the immersed boundary method (Peskin, 2002;Mittal and Iaccarino, 2005), or the lattice Boltzmann method (Chen and Doolen, 1998;Succi and Succi, 2001). However, the high resolution fluid flow and drag data associated with these simulations result in a high computational cost. Consequently, there have been relatively few documented studies within geomechanics using these fully resolved coupled models. Coupled, coarse-grained, unresolved simulations following the Anderson and Jackson (1967) formulation are more tractable (Kawaguchi et al., 1998;Kafui et al., 2002;Zeghal and El Shamy, 2004). Typically, these unresolved simulations use CFD to model the fluid phase coupled with DEM to model the solid phase, such that the fluid phase is discretised at a scale that is of the order of several particle diameters. However, such unresolved CFD-DEM cannot accurately capture local flow field effects and seepage induced drag, as quantities are averaged over a fluid cell comprising several particles. Furthermore, the drag coefficients commonly used in unresolved CFD-DEM in soil mechanics research do not account for the effect of polydispersity which are known to influence the fluid-particle interaction force (van der Hoef et al., 2005).
Pore Network Models (PNM) present an alternate numerical technique to model fluid flow through saturated sands. In these models, the nodes representing individual pores are connected by edges representing the narrow constrictions or pore throats. Conceptually, pore network models resolve the fluid flow at length-scale that is intermediate between the fully resolved and unresolved simulations approaches mentioned above. They have been employed in a wide range of fields including petroleum engineering (Silin and Patzek, 2006;Thompson et al., 2008), nuclear waste disposal (Xiong et al., 2015), ionic diffusion (Mohajeri et al., 2010) and carbon dioxide sequestration (Li et al., 2006). Within soil mechanics and geotechnical engineering, PNM have been employed to explore single-phase permeability (Bryant et al., 1993;van der Linden et al., 2018), unsaturated soil characteristics and fluid retention curves (Held and Celia, 2001;Ferraro et al., 2017) and filtration properties (Shire and O'Sullivan, 2017). The seepage induced drag force can also be obtained from PNM using the method developed in Chareyre et al. (2012).
While PNM are a computationally efficient approach to model fluid flow and fluid-particle interactions, they greatly simplify the geometry of the pore space and the complex physics associated with flow through a granular, porous medium (Miao et al., 2017). Network models are typically verified by comparing the macroscopic permeability with experimental data (Bryant et al., 1993) or fully resolved numerical simulations (van der Linden et al., 2018). The ability of PNM to capture local flow field effects has only been considered in a few studies. Helland et al. (2008) and Miao et al. (2017) compared local conductivity in PNM against finite element simulations using rock imaging data, while the lattice-Boltzmann method was employed by Sholokhova et al. (2009) for a sandstone sample. However, these studies were limited to considering local conductivity in isolated pore constrictions, and hence, did not capture the influence of neighbouring pores and constrictions on the local flow field behaviour. Low permeability rock samples were considered rather than granular "unconsolidated" sand samples. Moreover, PNM simplifications affect the ability to predict the seepage induced drag forces. The method proposed in Chareyre et al. (2012) was verified against finite element simulations for a simple cubic unit cell (comprising 8 particles) with an additional particle of variable size in the centre, while Catalano et al. (2014) extended the verification by comparing with the analytical solution for one-dimensional consolidation. Nevertheless, a more detailed validation of this method is required to establish whether PNM can accurately capture seepage induced drag. Therefore, this study aims to address the following key questions: 1. How accurately can a pore network model predict the local pressures at individual pores and the fluxes through individual constrictions across a broad range of packing fractions?
2. How sensitive are the results to the choice of local conductance models and the way the void space is partitioned?
3. How accurately can existing pore network models predict the seepage induced drag on individual particles?
The study compares the data generated from PNM with data from fully-resolved numerical simulations employing the Immersed Boundary Method (IBM) for a wide range of linearly graded and bimodal samples. The IBM simulations on entire samples capture the effect of neighbouring pores and constrictions, which is absent in prior studies. Further, the seepage induced drag determined from PNM is compared with analytical solutions for various ordered packings (including simple cubic, body centred cubic and face centred cubic arrays) provided in Zick and Homsy (1982). This analytical validation, in conjunction with the comparative study of drag forces obtained from IBM, provided much stronger support for the use of PNM to model fluid-particle systems.

Immersed Boundary Method
The equations of motion for an incompressible fluid in a multi-phase (fluid-particle) system are given by the Navier-Stokes equations: ∇ · u = 0 (1) where u is the fluid velocity, p is the pressure, f f p is the force to account for the presence of the local particle surface segments, ρ is the fluid density and ν is the kinematic viscosity. These equations express the mass and momentum conservation of the fluid. The solution to Eqs. 1-2 can be obtained using Computational Fluid Dynamics (CFD), a set of methods which numerically solves the Navier-Stokes equations by discretising the continuous fluid domain into a contiguous space-filling set of finite-volume fluid cells. The equations of motion for fluid cells form a system of non-linear equations, which are solved for the fluid velocity, u, and the pressure, p, at each fluid cell within the domain. Appropriate boundary conditions must be specified for fluid velocity and pressure at the domain boundaries, along with initial values for each fluid cell.
The Immersed Boundary Method (IBM) is a class of CFD methods in which the fluid equations are modified by the introduction of additional boundary conditions inside the fluid domain. Points representing the surfaces of solid objects (i.e. particle surfaces) are inserted within the fluid domain and no-slip, no-penetration boundary conditions are imposed by the addition of source terms representing the fluid-particle interaction force, f f p , in Eq. 2. This approach enables simulation of flow through the pore space in granular, porous media. In the IBM, the fluid mesh grid spacing is typically at least an order of magnitude smaller than the individual particles to enable the local flow field around the particles to be captured. The points defining the fluid cells are referred to as Eulerian points, as they remain fixed in space as fluid flows between them, while the mesh points representing particle surfaces are referred to as Lagrangian points, as they track the motion of the particles across the Eulerian fluid mesh. An illustration of the Lagrangian and Eulerian meshes is shown in Fig. 1. This study used MultiFlow (Denner and van Wachem, 2014), a fully-coupled fluid solver which solves the Navier-Stokes equations using the iterative biconjugate gradient stabilised method (van der Vorst, 1992) with a collocated variable arrangement. In this study, all IBM simulations were conducted with structured Cartesian grids. A multi-direct forcing type IBM is implemented in MultiFlow similar to that described by Luo et al. (2007). Assuming the particle properties are known at step N , at the next fluid time step, N + 1, the velocity of Lagrangian point l at coordinates X l i on the surface of particle i whose centroid is at a position x p i is given by where v i and ω i are the translational and rotational velocity vectors of particle i. In order to satisfy the no-slip condition at the next fluid time step, the force from Lagrangian point l to be applied to the neighbouring fluid at time step N is whereÛ N (X l i ) is an intermediate fluid velocity at X l i and time step N which satisfies Eq. 2 in the absence of the forcing due to the solid phase. The use of upper case symbols in Eq. 3 and 4 indicates variables defined at the Lagrangian points. Fluid velocities must be interpolated from the Eulerian mesh to the Lagrangian points such that where x m is the position of Eulerian cell m, δ h is the regularised Dirac delta function (DDF) of Roma et al. (1999), and ∆x, ∆y, and ∆z are the Eulerian mesh spacings in the x, y, and z directions. Correspondingly, forces calculated at the Lagrangian points must be distributed across all Eulerian cells within the support region of each point, which is three Eulerian cells wide for the DDF, such that where the sums are taken over the total number of immersed objects, N P , and the total number of Lagrangian points, N L , on each object and ∆V l i is the volume of the Lagrangian point. ∆V l i is related to the Lagrangian point area density and should match closely with the Eulerian cell volume for correct application of the DDF.
The requirement of no-slip between the particle surface and neighbouring fluid dictates that the fluid velocity at each Lagrangian point be equal to the rigid body velocity of the point as given by Eq. 3. Due to overlapping of the support regions, each fluid cell is influenced by multiple Lagrangian points and hence the values of u obtained at time step N + 1 will not satisfy the no-slip condition well. Multiple rounds of forcing are applied until the velocity field at time step N + 1 satisfies the no-slip condition to within a specified tolerance as measured by the difference between the target and interpolated fluid velocities at each Lagrangian point. The interpolation and spreading of variables also causes an enhancement of particle size from the perspective of the fluid, which leads to an overprediction of f f p . The radius retraction scheme (Breugem, 2012) is used to mitigate this effect.

Sample Generation
Samples were generated by producing a non-interacting cloud of spherical particles with a specified particle size distribution (PSD) and subsequently isotropically compressed to a target confining pressure in a cubic periodic cell. The particle coordinates following compression were used to construct the input for the IBM and PNM simulations.
This study considered linear graded samples with coefficient of uniformity, C u = 1.01, 1.20, 1.50, 2.00 ( Fig. 2) and bimodal samples with a particle size ratio of maximum to minimum diameter of χ = 1.99, 3.97 and various fines fractions (volume fraction of smaller particles) of ψ ≈ 0.10, 0.25, 0.35, 0.50. In all cases, the minimum particle diameter was taken as d min = 5.0 × 10 −4 m. A Monte-Carlo procedure was implemented in MultiFlow to generate the non-interacting particle assembly. The PSD was split into several size classes, and the number of particles and their diameters for each size class was calculated to match the target PSD. Particle diameters were distributed uniformly within each size class. Particles were placed within the domain in order of decreasing diameters, as the probability of placing larger particles is decreased if the space is already occupied by the smaller particles. Trial positions are only accepted if there were no overlap with other particles.
After the non-interacting particle assembly was generated, each sample was isotropically compressed (with no inter-particle friction) to a confining pressure of 100kPa using servo-controlled periodic boundaries and the MultiFlow DEM functionality. Snapshots were taken throughout the isotropic compression phase to obtain six configurations for each linear graded samples and two configurations for each bimodal sample. For the linear graded samples a wide range of packing fractions, φ = Vs V (V s is the total particle volume and V is the total sample volume), of interest to both chemical engineering and geotechnical engineering applications was considered. The least dense cases (φ ≈ 0.40) correspond to heavily laden particle suspensions found in conveying and materials processing applications, whilst the most dense cases (φ ≈ 0.65) correspond to percolating stress transmitting networks encountered in soil mechanics. For each sample, similar density values were considered to enable direct comparison. Table 1 lists the coefficient of uniformity (C u ), the ratio of maximum to minimum particle diameter (χ = d max /d min ), the number of particles (N p ), packing fraction (φ), and void ratio (e), of the linear graded samples considered in this study. The bimodal samples focused on stress percolating assemblies at two different densities in order to vary the PSD shape and explore the influence of fines content. Table 2 lists the particle size ratio (χ), fines fraction (ψ), number of particles (N p ), number of fine particles N F p , number of coarse particles N C p , packing fraction (φ) and void ratio (e). Following isotropic compression, the particle positions and diameters for each configuration were used as the initial particle positions for the IBM simulation. In the IBM simulations, the particles were fixed in position and a flow with low Reynolds number was allowed to develop to a steady state solution. Periodic boundaries were applied in the y and z directions only. In the x direction, the inlet boundary conditions (u x , u y , u z ) = (U, 0, 0) dp dx = 0 were applied at the lower end of the domain. The inlet velocity was set to U = 2 × 10 −4 m/s to simulate laminar flow (low Reynolds number, Re ≈ 0.4). The outlet boundary conditions were applied at the upper end of the domain. All the Lagrangian points on the particle surface must lie inside the domain. Placing Lagrangian points too close to the inlet will result in a non-physical solution as the no-slip condition imposed by these points will be in conflict with the constant velocity enforced by the inlet boundary condition.  Validation studies using ordered packings revealed that f f p was slightly over-predicted for particles at the inlet region and slightly under-predicted for those at the outlet region. In order to avoid these issues, the procedure illustrated in Fig. 3 was adopted, as explained below. Particles which crossed the DEM simulation boundaries along the x direction, are duplicated and translated to the opposite end of the domain. The domain is extended in the x direction, such that all Lagrangian points are separated by at least 2∆x from the inlet and outlet boundaries (where ∆x is the size of the fluid cell). When extracting the fluid-particle interaction forces, only those Lagrangian points which lie inside the original extents of the DEM simulation domain were considered. For a particle which is present at the inlet and also as a ghost particle duplicated at the outlet, part of the contribution to f f p is due to the flow field in the inlet region, and the contribution from the remainder of the Lagrangian points is due to the flow field in the outlet region. This procedure effectively mimics periodic boundary conditions along the x direction on u x , u y , u z and dp dx for the flow inside the original extents of the DEM simulation domain. The fact that f f p depends on the values of ∇p across the particles as opposed to absolute pressure, p, is a sufficient justification for the use of this approach. Tables 1-2 show that the assemblies considered in this study comprise very small number of particles. The IBM provides high resolution detail of the interstitial flow at great computational cost. To accurately capture details of the interstitial flow field, the smallest constrictions present must be spanned by a sufficient number of fluid cells. Here a resolution of dmin ∆x = 24 was adopted, where d min is the minimum particle diameter. As the range of particle sizes increases with increasing C u for linear graded samples or with increasing χ for bimodal samples, the overall domain size increases with more of the particle phase volume concentrated in a small number of large particles. This increase in domain size coupled with the mesh resolution criterion meant that the computational resources increased dramatically with C u and χ. This limited the sample sizes and particle gradations that could be considered in this study. A typical IBM case with approximately 600 particles embedded in an Eulerian mesh of 60 million fluid elements took approximately 16 hours to run on 400 processing cores. This imposed limits on the size and number of systems which could be considered.

Delaunay-based Network Generation
Pore Network Models (PNM) provide a simplified representation of flow in granular, porous media. The first stage in constructing a PNM is to partition the void space and identify the pores (PNM nodes) and pore throats or constrictions (PNM edges). Where assemblies of spheres are considered, a space tessellation method can be used, which identifies a set of unique, non-overlapping and spacefilling cells. The classical Delaunay Tessellation (DT) has been used in a number of prior studies of fluid flow in porous media (Mason, 1971;Bryant and Blunt, 1992;Thompson and Fogler, 1997). The DT produces a set of tetrahedral cells, where the vertices of the tetrahedron are the particle centroids. Each tetrahedron represents a node in the network model, while the faces of the tetrahedron are pore constrictions representing the edges in the network model. While the classical DT provides an adequate delineation of the pore space in monodisperse assemblies, it poses issues for assemblies of different sized spheres. Most notably, the faces of the tetrahedron may intersect adjacent particles (Chareyre et al., 2012). A weighted Delaunay Tessellation can overcome this issue by including consideration of the particle radius (r p ).
The weighted DT is applied to the samples generated for IBM simulations listed in Tables 1-2. This was conducted using the C++ package TetGen (Si, 2015), where the weighting considered was r 2 p . In order to avoid degenerate Delaunay cells that may develop on planar boundaries, a sub-sampling procedure was applied, taking into account the periodic domain employed in sample generation. In this sub-sampling process, the primary sample domain was replicated in all 26 adjacent directions and the DT was applied to this extended domain. Then, only those Delaunay cells completely within the primary domain were extracted and used in the network generation.
While the DT provides a simple and computationally efficient method to tessellate the pore space, Al-Raoush et al. (2003) suggested that it can overly subdivide the pore space and therefore not capture the existence of large pores. The DT also has a fixed network structure with all Delaunay cells comprising exactly four faces, and hence, four constrictions or edges in the network. To overcome these issues, Al-Raoush et al. (2003) suggested that tetrahedral Delaunay cells should be merged according to a merging criteria to form polyhedral pore cells, and termed this the Modified Delaunay Tessellation (MDT). The MDT results in a more physically representative tessellation of the pore space by capturing the large variation in pore sizes (that is, the existence of small and large pore bodies and constrictions) and allowing for variable edge connectivity in the network.
For the MDT considered here the merging criteria based on inscribed sphere (Al-Raoush et al., 2003;Reboul et al., 2008;Sufian et al., 2019) was used. An inscribed sphere of a Delaunay cell is the sphere which is internally tangent to the four particles forming the tetrahedron. A local optimisation procedure may be required if the calculated inscribed sphere is not contained entirely within the pore space. This is particularly prevalent in polydisperse assemblies. The local optimisation procedure outlined in Sufian et al. (2015) is employed in this study, which takes a combinatorics approach to determine the maximal inscribed sphere contained entirely within the pore space. Delaunay cells are merged if there respective inscribed spheres overlap.
The current study examines the sensitivity of PNM results to the way the void space is partitioned and so both DT and MDT data are used to construct the PNM. The geometric properties of the unit cells in the DT and MDT can be readily calculated analytically with the assumption of non-overlapping particles; an appropriate assumption for the simulations considered in this study. The geometric properties of interest are the pore volume, V p , the curved particle surface area within the unit cell, S p , and the pore area of the constriction, A p , which are used in the subsequent sections to define the network conductance model and calculate the fluid-particle interaction force. In the DT, S p is calculated from the spherical triangle formula, while V p is the determined by subtracting the volume of the particle sectors from the volume of the tetrahedral Delaunay cell, and A p is given by the area of the triangular face minus the area of the intersecting particle sectors (further details can be found in Sufian et al. (2015)). The polyhedral pores in the MDT are a union of non-overlapping and space-filling Delaunay cells. Hence, their geometric properties are given by the summation of corresponding Delaunay cells.

Local Pressure and Flow Field
Local fluid flow properties are obtained by simulating flow in the pore network using a Stokes flow algorithm. This is conducted using the open-source Python package, OpenPNM (Gostick et al., 2016). As in the IBM simulations, particles are fixed in position when conducting the flow simulation.
The Stokes flow algorithm ensures conservation of mass in the pore network, which at steady state conditions and in the absence of any source/sink terms (such as volume change in the pores) is given by: where q ij is the flux (flow rate) from node i to node j, g ij is the conductance of the ij edge, {p i , p j } are the nodal pressures, and the summation is over all n edges connected to node i. Models used to define the conductance of a network edge are detailed in Sec. 3.3. Application of Eq. 7 to all nodes in the network, taking the pressure at the inlet and outlet nodes as boundary conditions, forms a system of linear equations which can be solved using matrix-inversion algorithms to obtain the pressure at each node and consequently the flux along each edge. Dirichlet boundary conditions are imposed in the solution to Eq. 7, where pressure is prescribed at the inlet and outlet nodes. Inlet and outlet nodes are identified as those unit cells that have at least one face that is not shared with any other unit cells (which can only occur for those unit cells at the domain boundary) and are located in close proximity to the domain boundaries in the direction of flow. Note that an inlet flow velocity is considered for the boundary conditions in the IBM simulations. For consistency with the IBM simulations, the pressure at the DEM boundaries (indicated by the red dotted line in Fig. 3) are considered in the PNM to define the Dirichlet boundary conditions at the inlet and outlet nodes. A linear interpolation is performed to obtain the prescribed pressure at the centroid of the inlet and outlet nodes. The imposed pressure gradient in the PNM is equivalent and consistent with the imposed flow velocity (or flux) in the IBM, thereby enabling direct comparison between the two numerical techniques.

Conductance Models
The conductance of a network edge is defined by the geometry of the constrictions and/or pore bodies and the viscosity of the fluid. Two approaches are compared in this paper: (i) a hierarchical model based on pore-throat-pore series elements closely following van der Linden et al. (2018)

Pore-Throat-Pore Series Model
This model considers a simplified representation of the pore space as a collection of cylindrical elements connected in series and parallel. Flow through a cylindrical element is defined by the Hagen-Poiseuille equation: where q is the flux, r is the radius of the cylindrical element, l is the length of the cylinder, µ is the dynamic viscosity of the fluid, g is the conductance of the cylindrical element and ∆p is the pressure difference across the length of the element. The Hagen-Poiseuille equation is applicable for an incompressible fluid experiencing laminar flow, which is in line with fluid flow conditions in the IBM simulations. This conductance model represents the network edge element as a pore-throat-pore series, comprising two cylindrical pore elements connected in series by one or more cylindrical constriction element . The radius of the pore element, r p , is defined as the radius of an equivalent sphere of pore volume, V p . The length of the pore element, l p , is the length of an equivalent cylinder with volume, V p , and radius, r p , calculated as above.
In the DT, each pair of pores is connected by a single constriction, as there is only one shared triangular face between adjacent cells. However, multiple triangular faces may connect adjacent cells in the MDT. These can be regarded as multiple constrictions, where radii of the i th constriction is r t,i and the length is l t,i . The radius of the constriction elements, r t,i , is given by the radius of a circle whose area equals the pore constriction area, A p,i . The length of the throat element, l t,i , is given by the distance between the centroids of the two corresponding Delaunay cells to the shared triangular face. These multiple constrictions are connected in parallel, and can be expressed as a single equivalent element (Fig. 4b) by noting that flux in the single equivalent element is equal to the sum of flux in each constriction, q t = i q t,i and that the pressure drop is the same across all constrictions, ∆p t = ∆p t,i . If it is assumed that all constrictions connecting two pores have the same length defined by a weighted arithmetic mean, Figure 4: Schematic illustration of the conductance model: (a) two pore (node) elements connected by multiple constriction (edge) elements; (b) pore-throat-pore series element, where multiple throats are upscaled to an equivalent constriction; (c) network edge element formed by upscaling the pore-throatpore element.
then the radius of the single equivalent constriction element is given by: The pressure drop across a network edge element is the sum of pressure drops across each component in the corresponding pore-throat-pore series element, and the flux through the edge equals the porethroat-pore element flux. If it is assumed that a network edge element has a length, l c , defined as the centre-to-centre distance of the respective pore bodies (this differs from the centre-to-centre distance of corresponding Delaunay cells in the definition of l t,i for the MDT), then it can be shown that: where r c is the radius of the network edge element. Therefore the conductance of a network edge is given by:

Hydraulic Radius Model
While the pore-throat-pore series model in Sec. 3.3.1 used an hierarchical up-scaling technique to obtain an equivalent conductance for the network edge, the hydraulic radius model simply specifies the conductance of the network edge based on the geometry of the pores and constriction. In this model, the conductance is given by: where A t is the total pore constriction area (A t = i A t,i for the MDT), R h is the hydraulic radius, l c is the centre-to-centre distance of the respective pore bodies (as used in Sec. 3.3.1) and α is a nondimensional factor to account for the cross-sectional shape of constrictions. Following Chareyre et al. (2012), α = 0.5, which was deemed to be appropriate for assemblies of spheres. The hydraulic radius, R h , is defined as the ratio of volume to surface area of the two pore bodies corresponding to the shared constriction. Note that both conductance models (Eq. 11 and Eq. 12) utilise the geometry of pores and constrictions in defining the local conductivity and this is preferred over conductance models which exclusively consider constriction geometry (e.g. Bryant and Blunt (1992)) as the local flow field is influenced by the properties of both pores and constrictions. The conductance calculated from Eq. 12 is equivalent to the Hagen-Poiseuille form in Eq. 11 when r c = 2R h and α = 0.5. Therefore, the hydraulic radius approach with varying shape factors could be applied to each of the elements in the pore-throat-pore series and subsequently up-scaled to obtain an equivalent conductance. This is beyond the scope of this paper but demonstrates that the two conductance models are closely linked.
The conductance specified in Eq. 11 or Eq. 12 is supplied to OpenPNM and used in the mass conservation equation in Eq. 7 to determine the local pressure and flux in the pore network. Chareyre et al. (2012) proposed calculating the fluid-particle interaction force, F f p , by considering the geometric properties of the constrictions, along with the local pressure gradient obtained from the Stokes flow algorithm. F f p can be split into contribution of the pressure, F p f p , and viscous shear, F v f p , components (the subscript f p is omitted hereafter for clarity):

Fluid-Particle Interaction Force
where the superscript k indicates the force experienced by an individual particle. The pressure contribution to particle k from flow in element ij of the network is given by: where A k ij is the area of the sector of particle k intersecting the triangular face of constriction ij, {p i , p j } are the corresponding nodal pressures and n ij is the unit vector pointing from p i to p j . In a pore network representation, the pressure everywhere within a pore cell is assumed to be constant. In this study, the nodal pressures are defined at the cell centroid locations, which are located within the pore space.
The viscous shear contribution due to edge ij in the network is given by: where A f ij is the pore area in the ij constriction. The total viscous shear force is distributed to the particles intersecting the constriction by a ratio of their surface areas within the local domain: where S k ij is the surface area of the particle k within the local domain formed by nodes i and j, as illustrated in Fig. 5. The summation in the denominator to Eq. 16 is the total surface area within the local domain of the three particles forming the pore constriction. Figure 5: Schematic illustration of the terms required for the calculation of the fluid-particle interaction force. A k ij is the particle area intersecting the constriction, shown in red. A f ij is the pore or fluid area of the constriction, shown in blue. S k ij is the particle surface area within the union of the unit cells adjacent to the constriction, shown in orange. Note that the particle surface area for only one particle is shown for clarity.

Overall Flow Behaviour
A comparison of the overall flow behaviour is presented in Fig. 6, which shows the sample permeability obtained from the IBM and PNM. The PNM was generated using both the DT and MDT to explore the sensitivity to pore space partitioning. Both conductance models presented in Sec. 3.3 are also compared. The permeability, k, is obtained from the PNM by: where Q * is the summed fluxes at the inlet or outlet nodes, µ is the dynamic viscosity of the fluid, l is the domain length, A is the domain cross-sectional area and ∆p is the pressure difference across the domain length. Fig. 6 also plots the Kozeny-Carman estimate of permeability, defined as: where S v is the specific internal surface area to volume ratio obtained from the complete PSD by: where the summation is performed over all particles in the system. Fig. 6 demonstrates that the permeability obtained from the IBM simulations closely match the Kozeny-Carman expression. The permeability obtained from the PNM for the pore-throat-pore series conductance model (Sec. 3.3.1) is shown in Fig. 6a, while Fig. 6b shows the permeability obtained using the hydraulic radius model (Sec. 3.3.2). The hydraulic radius model better captures the overall flow response across the range of linear gradations and packing fractions considered in this study. The permeability obtained from both DT and MDT closely matches the IBM simulations, with slight overestimation at low packing fractions. In the pore-throat-pore series model, the permeability obtained from the PNM based on the MDT is in reasonable agreement with the value calculated in the IBM simulations, with permeability slightly underestimated. This suggests that partitioning the pore space more accurately using the MDT can result is a better estimate of permeability, irrespective of the conductance model. However, significant differences are noted for the permeability obtained from the DT when considering the pore-throat-pore series conductance model (see right-side figure in Fig. 6a), as the least dense assembly showed an order of magnitude error, with the error reducing with increasing density.

Correlation in Local Pressure and Flux
Local pressure and flux was obtained from the IBM simulations using the DT and MDT geometries that were input to the PNM, thereby allowing for direct comparison between the two numerical techniques. Pressures at the centroid of the DT cells were interpolated from the pressure fields obtained from the IBM simulations. A volume weighted sum of the pressures of the constituent DT cells was used to obtain each of the MDT cell pressures. These pressures calculated from the IBM simulations are compared against the nodal pressures obtained from the PNM. To calculate fluxes from the IBM data, each triangular face of the DT was discretised into a set of sub-triangles using a Delaunay triangulation such that the average area of the sub-triangles was approximately equal to the average area of the Eulerian cell faces. Fluid velocities were interpolated at the centres of each of these sub-triangles, allowing for the flux at each sub-triangle to be calculated. Fluxes through all sub-triangles were summed to obtain the total flux through each constriction in the DT. These DT fluxes were then summed to obtain the overall fluxes through the mutual faces between adjacent cells in the MDT. The IBM flux is compared against the flux through the corresponding edge in the PNM.     . For all packing fractions considered in this study, the local pressure field exhibited a near-perfect correlation between the IBM and PNM. This was evident for both the DT and MDT partitioning, along with both conductance models. This is further reflected in the Pearson correlation coefficient, ρ x , which is given by: where x is the property of interest (that is, pressure, flux or force), x PNM i is the value in the PNM, x IBM i is the value in the IBM, x PNM is the mean value of that property in the PNM, and x IBM is the mean value of that property in the IBM. The Pearson correlation coefficient has values ρ x ∈ [−1, 1], where ρ x = 1 indicates perfect correspondence between PNM and IBM values. In Table 3, the Pearson correlation coefficient for pressure (ρ DT p and ρ MDT p ) is above 0.99 for all linear graded samples, which confirms the visual observation of strong correlation. The correct identification of inlet and outlet nodes, along with the interpolated pressure at inlet and outlet nodes (as discussed in Sec. 3.2) was found to have a significant effect on the strong correlation in local pressure.
A reasonably strong correlation is maintained in the bimodal samples (Figs. 9-10), although slightly more scatter is observed. This is reflected in the Pearson correlation coefficients shown in Table 4, which remains greater than 0.97 for all bimodal samples. It is acknowledged that the bimodal samples with χ = 3.97 comprise very few coarse particles (N C p < 40 in Table 2), which may compromise the statistical significance of the results. Nevertheless, this computational limitation does not detract from the observation that the PNM can accurately capture the local pressure field irrespective of the method of tessellation or the specified conductance model.    The correlation in local flux is shown in Figs. 11-12 for linear graded samples using both conductance models, and in Figs. 13-14 for the bimodal samples. Unlike the near-perfect correlation in the local pressure field, the local flux correlation shows significant scatter.
When applying the pore-throat-pore series conductance model, Fig. 11 and Fig. 13 shows significant difference in the local flux correlation when the DT and MDT partitions are compared, with the MDT better capturing the local flux. In the linear graded samples, the difference between the MDT and DT predictions is most noticeable in the least dense samples (φ ≈ 0.40, first row in Fig. 11). The difference between the DT and MDT data reduces with increasing packing fraction, as evident for the cases with φ ≈ 0.65 (last row in Fig. 11). An increase in C u results in a wider scatter in the correlation plots for all packing fractions, which can be seen by comparing the first column in Fig. 11 for C u = 1.01 with the last column in Fig. 11 for C u = 2.00.
Similar trends are observed in the bimodal samples (Fig. 13), noting that the bimodal cases only consider stress-percolating assemblies (and have equivalent density to the last two rows in Fig. 11 for linear graded samples). As such, the difference in flux correlation for DT and MDT is less pronounced, but the MDT is consistently in better agreement with the IBM data. Furthermore, there is wider scatter for data points at large fluxes, particularly in the MDT. This is associated with non-physically large constrictions which have not been merged, as also recognised by Shire et al. (2013) in quantifying anisotropy of constrictions.
The difference in local flux correlation due to the DT and MDT partitions is less pronounced when considering the hydraulic radius conductance model ( Fig. 12 for linear graded samples and Fig. 14 for bimodal samples). The DT method exhibits slightly larger scatter, confirming that the MDT better captures local flux correlation. Moreover, the correlation appears to improve with increasing packing fraction. This can be visualised by comparing the first row (low packing fraction) with the last row (high packing fraction) in Fig. 12 for the linear graded samples. This differed from the observations for the pore-throat-pore series model in Fig. 11, which showed the opposing trend, suggesting that certain conductance models are more appropriate for certain ranges of packing fractions.
These visual observations of moderate correlation in local flux is confirmed by the Pearson correlation coefficient in Table 3 for local flux (ρ DT q and ρ MDT q ), which generally take a value between 0.6 and 0.8. Note that the clear visual differences in the DT and MDT for the pore-throat-pore series model ( Fig. 11 and Fig. 13) are not apparent when comparing the Pearson correlation coefficient, as ρ x quantifies how close the two quantities are to any linear relationship and not necessarily the linear relationship with unit gradient (as marked on Fig. 11), which corresponds to a perfect match between the IBM and PNM data. The Pearson correlation coefficients for the bimodal samples are shown in Table 4. In general, a wider scatter is observed with increasing particle size ratio, χ, and with decreasing fines content, ψ, but is of a similar magnitude to that seen in the linear graded samples. Figure 11: Correlation in normalised local flux (flow rate) in the PNM (q PNM ) and the IBM (q IBM ) for linear graded samples using the pore-throat-pore series conductance model in Sec. 3.3.1. Normalisation is conducted by q = q qave , where q is the flux through network edge element and q ave = U A x is the average sample flux with U = 2 × 10 −4 m/s is the boundary inlet velocity and A x is the cross-sectional area of the sample. Red markers indicate the data points for DT, while blue markers indicate data points for MDT. Figure 12: Correlation in normalised local flux (flow rate) in the PNM (q PNM ) and the IBM (q IBM ) for linear graded samples using the hydraulic radius conductance model in Sec. 3.3.2. Normalisation is conducted as per Fig. 11. Red markers indicate the data points for DT, while blue markers indicate data points for MDT. Figure 13: Correlation in normalised local flux (flow rate) in the PNM (q PNM ) and the IBM (q IBM ) for bimodal samples using the pore-throat-pore series conductance model in Sec. 3.3.1. Normalisation is conducted as per Fig. 11. Red markers indicate the data points for DT, while blue markers indicate data points for MDT. Figure 14: Correlation in normalised local flux (flow rate) in the PNM (q PNM ) and the IBM (q IBM ) for bimodal samples using the hydraulic radius conductance model in Sec. 3.3.2. Normalisation is conducted as per Fig. 11. Red markers indicate the data points for DT, while blue markers indicate data points for MDT.

Shortest Paths and Streamline Profiles
Flow in a granular, porous medium is dominated by preferential flow channels. To evaluate the ability of PNM to capture these flow channels through the samples predicted by the IBM simulations, the shortest paths between pairs of inlet and outlet nodes were determined using Dijkstra's shortest path algorithm as implemented in the C++ software library igraph (Csardi and Nepusz, 2006). A weighted, directed graph was constructed based on the PNM. Directionality was assigned by considering the fluxes, q ij , along the edges of the graph, where edges with q ij < 0 were reversed to ensure that all edges experienced a positive flux. A weight, w ij was assigned to each edge and was set to be inversely proportional to the flux, q ij , in the corresponding constriction. The shortest path is the path which minimises the sum of the weights on the edges connecting an inlet and outlet node. This effectively captures the paths along which tracer particles would be transmitted most quickly. The shortest path from every inlet node to every outlet node was determined. Fig. 15 shows the ten shortest paths in the entire assembly based on the weighted, directed graph obtained from IBM and PNM fluxes. Despite the differences in local fluxes observed in Figs. 11-13, it is clear that the shortest paths predicted by the IBM and PNM are in close agreement.
The shortest paths can be described as ordered sequence of indices of the nodes traversed along the path. To quantify the similarity of the shortest paths demonstrated visually in Fig. 15, a measure of the similarity, S, of pairs of paths between nodes indices a and b was defined as: An example is provided to explain the above path similarity measure. For a path starting at inlet node 0 and ending at outlet node 5, the set of nodal indices the path passes through in the IBM results may be 0 P 5 IBM = {0, 1, 2, 3, 5}, and in the PNM results may be 0 P 5 PNM = {0, 1, 2, 4, 5}. In this case, the similarity of the pair of paths from the IBM and PNM results between nodes 0 and 5 has a value S(0, 5) = 4/5. This is a simplified example; flow paths in the present study typically traverse between 10 and 20 nodes.
For each inlet node, the shortest N SP paths are obtained from the IBM and PNM results. The N SP pairs of paths with common outlet nodes were identified and a mean value of S was calculated for each inlet node: where b j is the index of the outlet node of the j th path of the N SP paths from the input node with index a which have common output nodes in the IBM and PNM data. The probability density functions and cumulative probability distributions of S for the densest case considered for the C u = 1.01, 1.20, 1.50, and 2.00 cases calculated with N SP = 90 and N SP = 10 are shown in Fig. 16 and Fig. 17, respectively. The distributions of S shows that around half of the paths through the samples have in excess of 80% of the nodes in common between the IBM and PNM results for all values of C u considered. This provides strong quantitative confirmation that the flow paths predicted by the PNM match very closely with those predicted by the IBM, supporting what is visually observed in Fig. 15. Further, Fig. 17 shows that if only the top ten shortest paths are considered for each input node, then a better agreement is found, with the peak of the probability density function at S = 1 for all C u . In addition to shortest paths, streamlines provide further insight into dominant flow channels. A streamline is a parametric curve which is instantaneously tangent to a vector field at all points. Streamlines tracing the paths of massless particles through the interstitial space of the samples were extracted from the highly resolved flow fields obtained from the IBM simulations. Seed points specifying the starting location of the streamlines were uniformly distributed over the inlet boundary. At steady state flow conditions a streamline is fully defined by its seed location.
The role of pore and constriction geometry in pressure loss at the micro-scale was assessed by combining the extracted streamlines with the pore geometry and topology specified by the MDT in the PNM. Streamlines are parametrised in terms of time, t. Starting from the seed location at time, t = 0, as the tracer particle travels through the voids, it is enclosed within a single tetrahedral cell of the DT, and  Fig. 19. The colour of the streamline indicates the proximity to the labelled constrictions, as defined by C p (t) in Eq. 22.
hence a single pore cell of the MDT, at all times. The distance, C d (t), along the streamline to the closest constriction plane was determined for every point comprising each streamline. A constriction proximity measure, C p (t), was defined at each point along the streamlines: where α is a characteristic length with value α = 0.1 mm. This value was selected to obtain stiffly peaked Gaussians at the locations of the constrictions. In Fig. 18 a streamline extracted from the C u = 1.01 and φ = 0.651 case is shown. The traversed pore cells of the MDT are shown beneath the streamline, coloured by their cell indices to clearly delineate the boundaries between adjacent pore cells. The streamline is coloured according to the value of C p , such that points where the streamline crosses the constrictions are shown in red. Each constriction has been labelled. The velocity profile through each constriction is shown in Fig. 19. The pressure, p, and magnitude of the pressure gradient, |∇p|, along with the value of C p at each point along the streamline in Fig. 18 are shown in Fig. 20. Peaks identified in the C p data are marked across the three axes with black lines. Peaks in |∇p| are clearly observed to occur either exactly at, or very close to, a peak in C p , that is, at the constriction planes. The first peak and the two last peaks in |∇p| do not correspond to peaks in C p because these peaks occur at the boundary layers of the IBM data that are outside the domain of the MDT. In Fig. 21 pressure drop data is shown for the C u = 1.50, φ = 0.621 case which exhibits the same behaviour seen in Fig. 20. This feature is exhibited in all other cases (not presented here). Indraratna et al. (2012) suggested that the constriction size distribution rather than the particle size distribution is the dominant factor in determining the macro-scale pressure drop. Taylor et al. (2017) conducted CFD simulations with soil samples re-created from X-ray microcomputed tomography images of real sands and found that the majority of pressure loss occurred at the constrictions. The data presented here support these earlier findings and highlights the importance of the constriction geometry in pore scale pressure loss. These results also demonstrate the ability of the Modified Delaunay Tessellation based on the inscribed sphere procedure described in Sufian et al. (2015) to identify the constrictions which govern the permeability of porous media.   Figure 21: Pressure profile (a) and pressure gradient profile (b) against arc distance, d s , along the path of a stream tracer particle for the C u = 1.50, φ = 0.621 case. In the lower panel (c), C p , indicates the proximity to the nearest constriction.

Validation and Correlation for Fluid-Particle Interaction Force
A rigorous validation of the method proposed by Chareyre et al. (2012) as outlined in Sec. 3.4 is presented here by comparing against the analytical solution in Zick and Homsy (1982) for Stokes flow in various ordered packing configurations. The analytical solution for the (magnitude of the) fluid-particle interaction force is given by: where K represents a scalar multiple of the Stokes drag force experienced by a single particle, and U is the fluid velocity, which can be calculated from the pressure gradient across the length of the sample by: Zick and Homsy (1982) tabulated values for the factor, K, for various ordered configurations and packing fractions. In this study, validation is performed for simple cubic (SC), body centred cubic (BCC) and face centred cubic (FCC) configurations, with five packing fractions considered for each configuration (Fig. 22). In all cases, the particle radius was r p = 0.001m and assemblies of different densities are formed by increasing the lattice spacing. The ordered assemblies were generated by replicating the crystal lattice for each ordered configuration 4 × 4 × 4 in each direction. The number of particles (N p ) in each assembly can be found in Tables 5-7. Table 5: Comparison of fluid-particle interaction force in a simple cubic assembly from PNM simulation (K DT and K MDT ) and analytical solution (K SC ) using the pore-throat-pore (PTP) and hydraulic radius (HR) conductance model  191 39.205 44.014 48.623 42.191 42.191 42.191 44.625 47.213 49.745 42.191 42.191 42.191 Table 6: Comparison of fluid-particle interaction force in a body centred cubic assembly from PNM simulation (K DT and K MDT ) and analytical solution (K BCC ) using the pore-throat-pore (PTP) and hydraulic radius (HR) conductance model  Table 7: Comparison of fluid-particle interaction force in a face centred cubic assembly from PNM simulation (K DT and K MDT ) and analytical solution (K FCC ) using the pore-throat-pore (PTP) and hydraulic radius (HR) conductance model A pore network based on the DT and MDT was generated following exactly the same procedure as the linear graded and bimodal samples using both conductance models detailed in Sec. 3.3. The number of unit cells in the DT (N DT c ) and MDT (N MDT c ) is listed in Tables 5-7. In the SC cases, a significant amount of merging has occurred, but this remains consistent for all densities considered. In the BCC cases, no merging has taken place for all assemblies except for the least dense case, φ = 0.400, (a) Pore-throat-pore series model (b) Hydraulic radius model Figure 22: Validation of network-based calculation of the normalised fluid-particle interaction force, for ordered assemblies, including simple cubic (SC), body centred cubic (BCC) and face centred cubic (FCC). F f p is normalised by the Stokes drag force, which as per Eq. 23, is given by the scalar factor, K. Note that some of the DT markers are exactly in the same position as the MDT markers, and hence, cannot be seen in the figure. which due to the larger lattice spacing would allow for inscribed spheres to intersect. Some merging is observed in the FCC cases, but unlike the SC and BCC cases, the assemblies at φ = 0.460 and φ = 0.700 exhibit different numbers of unit cells in the DT. This is due to boundary artefacts with the Delaunay triangulation algorithm, where additional tetrahedral cells are detected at the boundaries for the FCC configuration. In this validation study, these boundary artefacts have not been removed, to ensure that the procedure followed in the validation study is exactly the same as that undertaken for the linear graded and bimodal samples, where this issue was not observed.
Flow is simulated in these ordered networks by applying an inlet pressure of 1.00Pa and an outlet pressure of 0.01Pa. Fig. 22 compares the numerical results against the analytical solution for various ordered configurations and packing fractions, where the factor, K, is obtained from a cubic spline interpolation of the tabulated values in Zick and Homsy (1982) and compared against the numerical value obtained from the PNM, K PNM , which is equal to the fluid-particle interaction force normalised by the equivalent Stokes force. The comparison is performed for only internal particles, as particles adjacent to the domain boundary have an incomplete force computation. The fluid-particle interaction force obtained from the MDT is more accurate compared to the DT. The only exception is the BCC case with φ = 0.400, which is due to the increased degree of merging not evident in other BCC cases. Also, the increase number of unit cells in the FCC cases at φ = 0.460 and φ = 0.700 only resulted in a very slight discrepancy, providing further justification that the boundary artefacts had minimal influence.
While Fig. 22 shows the average fluid-particle interaction force over all particles, further details are provided in Tables 5-7. In the SC cases, a range of forces is obtained for DT, as indicated by the different minimum (K DT min ) and maximum (K DT max ) values. For the pore-throat-pore conductance model, the fluidparticle interaction force was slightly underestimated (K DT mean < K SC ), while it was slightly overestimated for the hydraulic radius model when applying the DT. The same value is obtained for all particles when applying the MDT (as K MDT min = K MDT max = K SC , where K SC is the analytical solution based on Zick and Homsy (1982)). The same force value for all particles is also obtained in the BCC cases (with the exception of φ = 0.400 for the reasons discussed above). The FCC cases exhibit similar response with a range of forces obtained for DT and similar forces in MDT (there is only a slight difference in K MDT min and K MDT max for φ = 0.560, 0.640, 0.740). As per the discussion above, differences are noted for the FCC cases for φ = 0.460, 0.700. Therefore, the above validation provides strong support for the use of the method developed in Chareyre et al. (2012) to calculate the fluid-particle interaction force. While the DT provides a reasonable estimate, the use the MDT results in a more accurate prediction of the force. This provides confidence in applying the method to linear graded and bimodal samples.
Comparison of normalised fluid-particle interaction force, F k f p , obtained from the PNM and IBM for the linear graded samples is shown in Fig. 23-24 for both conductance models, along with DT and MDT partitions. This demonstrates that the PNM can reasonably capture the fluid-particle interaction force with most data points within a ±30% error bound, where the force obtained from the PNM is generally lower than that observed in the IBM. In Figs. 23-24, all cases are plotted on a single figure to show consistent trends between the samples, and Table 3 provides the Pearson correlation coefficient data. For the pore-throat-pore series conductance model, the MDT provides a more accurate agreement between the IBM and PNM, as reflected by ρ MDT F ≥ ρ DT F in Tables 3 with the correlation improving with increasing packing fraction. This is in line with observations for local pressure and flux correlations with this conductance model. In the hydraulic radius conductance model, both DT and MDT show similar degrees of correlation, with DT exhibiting slightly higher Pearson correlation values. Table 3 appears to suggest that the correlation improves within increasing C u . However, this is an artefact of the definition of the Pearson correlation coefficient in Eq. 19. At small C u , there is only a small particle size ratio and hence, only a small range in F k f p . Consequently, small differences can have a significant affect on the Pearson correlation coefficient. This should be considered when interpreting the correlation shown in Tables 3-4. This could be resolved by performing simulations (PNM and IBM) at various pressure gradients (using the same particle configuration), which would in turn result in a fluidparticle interaction force over wider range. This was not possible within the computational limits of this study, but it would be anticipated that the Pearson correlation coefficient would improved substantially for lower C u .
The calculation of the fluid-particle interaction force using the method proposed by Chareyre et al. (2012) relies on the pressure gradient across a constriction, along with the geometric properties of the constriction. In order to further test the simplifications and assumptions associated with the method, the normalised fluid-particle interaction force is recalculated by directly applying the nodal pressures obtained from the IBM simulations into the equations in Sec. 3.4 (instead of the calculated pressures obtained from the mass continuity equations in Eq. 7). The correlation is shown in Fig. 25 and the corresponding Pearson correlations values are included in Table 3 (see 'IBM' columns in the 'Force'  table). Fig. 25 shows a strong correlation, particularly with the MDT partition, with data points well within the ±30% error bound, which is confirmed by the Pearson correlation coefficients with value Normalisation is conducted by F = F F stokes , where F is the interaction force on a particle and F stokes = 3πµdU is the Stokes drag force on an individual particle (µ is the dynamic viscosity of the fluid, d is the particle diameter and U is the inlet velocity). The dashed-dotted line represents 30% deviation from perfect correspondence between IBM and PNM.  bound shown in Figs. 26-27 for both tessellation methods (DT and MDT) and conductance models. The large scatter is most prevalent in the χ = 3.97 samples. This is not reflected in the Pearson correlation coefficient in Table 4. In bimodal samples, the fine particles experience lower force, while the large particles experience higher force. This results in a very large variance in force, and consequently only very large errors will affect the correlation. Therefore, these values should be interpreted carefully. Following Fig. 25, the fluid-particle interaction force is recalculated using the IBM nodal pressures in the bimodal samples. The force correlation is shown in Fig. 28. Quite significant fluctuations are observed for the cases with χ = 3.97, which is a result of the relatively small size (particularly with the insufficient number of coarse particles) of these samples. These fluctuations are shown to alert readers of the need to have a representative element volume, which can increase significantly with increasing particle size ratio. By isolating the results for χ = 1.99, it is apparent that the method developed in Chareyre et al. (2012) can reasonably handle bimodal (and more generally gap-graded) samples. Nevertheless, the scatter observed in Fig. 28 for the bimodal cases is significantly more prominent that the minimal fluctuations observed in the linear graded samples (Fig. 25).

Concluding Remarks
This study addressed the ability of Pore Network Models (PNM) to capture the local pressure and flow field in granular, porous media, and the seepage induced drag (also referred to as the fluid-particle interaction force). Fully-resolved flow simulations were conducted using the Immersed Boundary Method (IBM) to obtain local pressure and flow fields at a sub-particle spatial resolution, along with the seepage induced forces acting on individual particles. The IBM simulations were conducted on linear graded samples with the coefficient of uniformity, C u , ranging from 1.01 to 2.00, and a wide range of packing fractions were considered for each gradation. In addition, bimodal samples were considered with a particle size ratio of 1.99 and 3.97, and at various fines content, ranging from 0.10 to 0.50.
Using these linear graded and bimodal samples, PNM was generated using Delaunay-based methods. The weighted Delaunay tessellation (DT) was considered, which accounted for the different sized particles by using a weighting factor based on the particle radius. This weighted DT was essential in properly delineating the pore space in the C u = 2.00 linear graded samples along with all bimodal samples, and hence, is the recommended approach for polydisperse sphere assemblies. In addition, the modified  Delaunay tessellation (MDT) was considered, which merged tetrahedral cells according an inscribed sphere criterion to form polyhedral unit cells. The MDT better captures variation in pore and constriction sizes, along with removing the connectivity constraint in the DT. Fluid flow in the network model was simulated using a Stokes flow algorithm, which ensured conservation of mass in the pore network. Two local conductance models were compared: (i) pore-throat-pore series element; and (ii) equivalent hydraulic radius approach. These conductance models considered the geometry of pores and constrictions in assigning a local conductivity to the network edge elements. The pressure at each node in the network, along with the flux (or flow rate) along each edge of the network was determined using the Stokes flow algorithm, with boundary conditions matching the fully-resolved IBM simulations. Comparison of the sample permeability showed that the PNM based on the hydraulic radius conductance model was able to closely match the permeability obtained from the IBM for both DT and MDT partitions. In the pore-throat-pore series conductance model, the MDT accurately captured the sample permeability, while the DT showed significant differences. The local pressure and flow field were also compared. Across the wide range of packing fractions considered in both linear-graded and bimodal samples, the local pressure field could be accurately captured by the PNM using either conductance models and for both the DT or MDT. However, the local flux showed more sensitivity to the choice of the conductance model and the method of partitioning the pore space (i.e. DT or MDT). In the pore-throat-pore series conductance model, local flux was reasonably captured using the MDT, while the DT showed significant differences with the IBM. In the hydraulic radius conductance model, DT and MDT methods exhibited similar behaviour (with the MDT being slightly more accurate).
Preferential flow channels were compared in the IBM and PNM by obtaining the shortest paths from the inlet to outlet nodes. A high degree of correspondence was observed between the IBM and PNM shortest paths, despite the differences in local fluxes. As flow characteristics of a granular, porous medium are dominated by the shortest paths, as they carry the majority of the flux of fluid through the sample, these results provide strong evidence supporting the use of PNM. Furthermore, streamline profiles demonstrated that pressure drops in the granular, porous medium coincided with the pore constrictions identified in the MDT, providing further evidence that the MDT provides a better representation of the pore space and should be considered when generating pore networks.
In addition, the fluid-particle interaction force was obtained from the PNM using the methodology proposed in Chareyre et al. (2012). This method was rigorously validated by comparing against the analytical solutions in Zick and Homsy (1982) for simple cubic, body centred cubic and face centred cubic array of spheres, which showed that the MDT was essential in accurately capturing the fluidparticle interaction force, although the DT provided a reasonable estimate. Following this validation, the force obtained from the PNM and IBM were compared for the linear graded and bimodal samples. The PNM provided a reasonable estimate of the interaction force for linear graded samples with most particles having an error bound of ±30%. However, the comparison was less promising for the bimodal samples, which exhibited a wider scatter. Interestingly, calculation of the fluid-particle interaction force was relatively insensitive to the network generation method with the MDT being slightly more accurate than the DT.
The differences in the local flux field could be attributed to the conductance model implemented, which may be adjusted to provide a better correlation. For example, different cross-sectional shapes could be considered in the hydraulic radius model, including a popular approach known as the circle-trianglesquare (CTS) method (Patzek and Silin, 2001), along with the consideration of adjustable parameters in the pore-throat-pore series conductance model as suggested by Sholokhova et al. (2009). Nevertheless, the choice of this conductance model is quite arbitrary, and this remains an unresolved feature of modelling pore networks. The findings in this study suggests that conductance models calibrated to match the overall permeability will reasonably capture the local pressure and flow field, as is the case for the PNM based on the MDT, but further research is required to support this statement. An alternate approach is to use machine learning techniques to obtain a conductance model without shape simplifications, as suggested by Miao et al. (2017), and this is a promising direction to resolve the issues surrounding the choice of the conductance model.
Differences in the calculated fluid-particle interaction force with the IBM results (particularly in the bimodal samples) may be attributed to the pore space discretisation, choice of the conductance model, or the method proposed in Chareyre et al. (2012). The force computation was relatively insensitive to the choice of the partitioning method (i.e. DT or MDT) or the choice of the conductance models considered in this study. The methodology in Chareyre et al. (2012) was further tested by directly applying the IBM nodal pressures in the force calculation, which exhibited excellent correlation in the linear graded samples. However, the bimodal samples exhibited some fluctuations despite the fact that the local pressure field was well captured. Therefore, this study provides strong support for the method developed in Chareyre et al. (2012) which is capable of capturing the fluid-particle interaction force, but further research is required to improve the accuracy of the method (particularly in bimodal or gap-graded particle assemblies) by further investigating discretisation methods, local conductance models and the underlying assumptions and simplifications in the force computation.