2.1. Grain-size distributions and sediment samples

In the laboratory, we collect data on grain-size distribution, porosity and hydraulic conductivity for both idealised and Earth materials: For our idealised materials, we create mixtures of glass shot and sand that we then compare them with natural subglacial tills.

To have a benchmark for our modelling, we aim to produce bi-modal distributions of two discrete grain sizes separated by approximately one order of magnitude. The actual separation we obtain is limited by the resolution of the mesh sizes bracketing each retained sample. We use spheroidal glass shot passed through a 180 ${\unicode[Arial]{x03BC}}$ m mesh (No. 80 sieve) that was retained on a 150 ${\unicode[Arial]{x03BC}}$ m mesh (No. 100 sieve). We mix this material with well-rounded quartz beach sand, sourced from the Bayfield Peninsula, Wisconsin, passed through a 2 mm mesh (No. 10 sieve) and retained on a 1 mm mesh (No. 18 sieve). Ultimately, we produce two approximately normal populations with mean diameters centred at 160 ${\unicode[Arial]{x03BC}}$ m $(\sigma _{F} = 29.7$ ${\unicode[Arial]{x03BC}}$ m $)$ and 1500 ${\unicode[Arial]{x03BC}}$ m $(\sigma _{C} = 320$ ${\unicode[Arial]{x03BC}}$ m $)$ for the shot and sand, respectively. We then mix these two contrasting materials together in three different proportions of volume – approximately 1 : 10, 1 : 1 and 10 : 1 shot to sand.

We measure the grain diameters for each idealised material using laser diffraction analysis (Malvern Mastersizer 2000E, range 0.01–2000 ${\unicode[Arial]{x03BC}}$ m). Before analysis, we wash and dry the samples to remove fines and prevent aggregation, then disperse them in deionised water. For each material, we perform three replicate runs, reporting our results as the volume-based median diameter ( $d_{50}$ ). For coarse sand, the particle sizes approach the upper detection limit of the instrument, making laser diffraction unreliable. Instead, we estimate grain diameters through image analysis. We photograph a representative subsample ( $n \approx 200$ grains) on a lightboard under uniform illumination, with scale calibration provided by a stage micrometer. We then process the images using ImageJ software to calculate equivalent circular diameters, accounting for the irregularity of grain shapes.

Figure 2.

Image of our mixture of two idealised synthetic materials alongside component grain-size distributions (top row) and image of subglacial till alongside grain-size distribution (bottom row). (a) Image of our mixture of glass shot and sand grains. (b) Grain-size distribution for fine glass shot. (c) Grain-size distribution for coarse sand grains (bottom). Relative frequencies are represented as bars, with best fit normal distributions as opaque dotted lines. (d) Image of subglacial till sampled from a deglaciated paleo deposit of the Bindschadler Ice Stream in West Antarctica. (e) Corresponding grain-size distribution for Bindschadler subglacial till.

Figure 2. Long description

Panels a-c show an idealized mixture, a close-up of fine glass beads, and histograms with normal fits for fine beads near 162 micrometers and coarse sand near 1,460 micrometers. Panels d-e show an image of Bindschadler Till and its broad, bimodal grain-size distribution. The inset enlarges the fine mode on a logarithmic grain-size axis.

We also conduct hydraulic conductivity experiments on a suite of natural tills from Antarctica. Ten of these tills were collected in the field and are now archived at The Polar Rock Repository (PRR) et al. (2015), along with one subglacial till sampled from a paleo–grounding zone of the Bindschadler Ice Stream (McGlannan et al. Reference McGlannan, Bart, Chow and DeCesare2017). The PRR samples represent a range of terrains and lithologies with further details for each sample available in the PRR database under the corresponding sample names. Although the precise genesis and depositional history of each sample are uncertain, they share key characteristics of basal tills, most notably a wide variance in grain sizes. For analysis, we remove clasts larger than 2 mm and then measure grain-size distributions of the processed sediments by laser diffraction using a Malvern Mastersizer, following standard protocols.

In figure 2, we compare the idealised mixture with the real Earth material both in appearance and in grain-size distribution. An example mixture of idealised materials is shown in figure 2(a), with an inset highlighting the fine glass shot. Histograms of grain sizes for the fine glass shot and coarse sand are shown in figures 2(b) and 2(c), respectively, along with fitted normal distributions obtained using maximum likelihood estimation. The fine glass shot histogram is coloured light blue, distinguished from the coarse sand histogram coloured in black. Both histograms and Gaussian fits are tightly centred on their respective modes, with normalised standard deviations (SDs) of $\sigma _{F} / \mu _{F} = 0.18$ for the fine mode and $\sigma _{C} / \mu _{C} = 0.21$ for the coarse mode. The two modes are separated in scale by roughly one order of magnitude, with $\mu _{F} / \mu _{C} = 0.11$ .

We show an example scanning electron microscopy (SEM) image of the subglacial till from Bindschadler Ice Stream (McGlannan et al. Reference McGlannan, Bart, Chow and DeCesare2017) in figure 2(d). The large image highlights some of the coarse sand grains and the inset highlights the extremely small and heterogeneous clay grains. We plot the corresponding probability distribution function of the grain sizes in figure 2(e), with an inset that represents the magnitude of the fine clay mode. Like the distinct colours used in figures 2(b) and 2(c), light blue is used in the inset of figure 2(e) to distinguish the fine mode from the coarse mode. The magnitudes and widths of the fine mode and coarse mode differ greatly for the till, making it challenging to assess the relative amount of fine grains and coarse grains visually. By directly computing the areas under the curve for the blue region and white regions, we find that the fraction of fine grains ( $0.68$ ) is of comparable magnitude to the fraction of coarse grains ( $0.32$ ).

When comparing the grain-size distributions in figures 2(b,c) and 2(e), the two modes of the idealised mixture are very different from the two modes of the subglacial till. In the idealised mixture, both the fine and coarse modes exhibit little variance relative to their mean grain size, the modes are separated by roughly an order of magnitude and the distributions for each mode do not interfere with each other at all. However, the subglacial till exhibits less clear separation: the fine and coarse modes see wide variances relative to their mean grain size, the modes are not very separated and the distributions for each mode clearly interfere with each other. In figure 3, we plot grain-size distributions of several PRR tills with an inset for the large peak for the mode at around 0.1–10 ${\unicode[Arial]{x03BC}}$ m. Subglacial tills have varying degrees of overlap among the two primary modes of grain sizes, and sometimes have additional modes (such as PRR 23842) or only a single distinguishable mode (such as PRR 23943).

2.2. Measuring hydraulic conductivity

We measure hydraulic conductivity in a four-inch diameter Gilson $^{\text{TM}}$ rigid-wall compaction-mould permeameter using the constant head method, following ASTM D5856-15 (ASTM International 2015). The only exception is Bindschadler Till, which has high clay content and low conductivity, and could not be tested in the rigid-wall device. Instead, we use a six-inch diameter flexible-wall permeameter with the falling head method, following ASTM D5084, to measure conductivity at three different porosities. For the idealised mixtures and the Bindschadler Till, we measure conductivity at multiple porosities, whereas each PRR till was tested at a single porosity. We achieve variable porosities by compacting samples in the appropriate mould to a target density. We estimate the resultant porosities from water content, mass and density using procedures outlined in the respective ASTM standards. Reported values of hydraulic conductivity from the constant head tests represent the mean of six measurements collected during a single experiment, while flexible-wall permeameter results represent averages of multiple measurements collected over 1–2 weeks of run time.

Because the material samples contain clay, the influence of clay swelling during fluid injection for hydraulic-conductivity measurements cannot be neglected (Scheuerman & Bergersen Reference Scheuerman and Bergersen1990). Tap water was used as the probing fluid, which should have reduced the likelihood of swelling relative to very low ionic-strength water. However, the ionic composition was not independently measured or controlled. Because several samples contain appreciable montmorillonite, fluid–clay interactions may still have influenced the measured hydraulic conductivities by altering pore-structure during measurement. Accordingly, the reported values are interpreted as effective conductivities under the experimental conditions employed, rather than intrinsic values independent of pore-fluid chemistry.

Figure 3.

Grain-size distribution for several subglacial tills. Insets depict the fine mode at approximately 0.1–10 ${\unicode[Arial]{x03BC}}$ m.

We report corresponding permeabilities because permeability removes the direct dependence of hydraulic conductivity on fluid density and viscosity, although it does not remove changes to the medium caused by fluid–solid interactions. Because our laboratory measurements are conducted using a constant head method with water at room temperature, we compute an effective permeability $k$ from the measured hydraulic conductivity $K$ using the relation:

(2.1) \begin{align} k = \frac {\mu }{\rho g}K , \end{align}

where we use $\mu = 1.0\times 10^{-3}$ Pa s as the dynamic viscosity of water at $20\,^\circ\textrm{C}$ , $g = 9.8$ m s $^{-2}$ as the acceleration due to gravity and $\rho = 1.0 \times 10^3$ kg m $^{-3}$ as the density of water at $20\,^\circ\textrm{C}$ (Kestin et al. Reference Kestin, Sokolov and Wakeham1978). We list these measured hydraulic conductivities and effective permeabilities for the subglacial tills from The Polar Rock Repository (PRR) et al. (2015) in table 1. While these subglacial tills have low permeabilities relative to other granular materials, tills have high permeabilities relative to many rocks that have been the object of many observations (Mahesar et al. Reference Mahesar, Memon, Memon and Tunio2017; Scott et al. Reference Scott, Lévy, Covell, Franzson, Gibert, Valfells, Newson, Frolova, Júlíusson and Guđjónsdóttir2023; Vali & Haji Zadeh Reference Vali and Haji Zadeh2024). In particular, these tills have permeabilities above the 1 mD, or $9.8 \times 10^{-16}$ m $^2$ , limit where the choice of interstitial fluid used for measurement can lead to differences of orders of magnitude (Mahesar et al. Reference Mahesar, Memon, Memon and Tunio2017).

Table 1.

Summary of laboratory measurements including grain-size distributions, as well as hydraulic conductivities and corresponding bulk permeabilities. Grain size data for idealised materials as well as more complete data for different tills from the Polar Rock Repository.

3.1. Conceptualising the dual porosity space

Some properties of granular materials, such as grain-size distribution, can be measured directly. Other properties, such as porosity, represent material averages and can only be measured in bulk. Bulk porosities capture the total amount of open space for water to travel through in a sample of material, but do not reveal how the open space is spatially distributed. At one extreme, a material with a given porosity could have a small number of large pores on a coarse length scale. At the other extreme, a material with the same porosity could have a large number of small pores at a fine length scale.

To address the variability of pore sizes in our model, we simplify the distribution of pore sizes into a dual-porosity system with pores of coarse-scale size and pores of fine-scale size. We then use the two modes of the grain-size distribution to estimate how much of the material likely consists of sparse coarse-scale pores or dense fine-scale pores. Using these data, we reduce some of the rich variability expected in the pore geometry to a deterministic system for our study.

Figure 4 shows our conceptualisation of the multi-scale pore geometry. To understand the dynamics of a polydisperse sediment layer, $\varOmega '''$ , with thickness $l'''$ , we study a unit cell $\varOmega ''$ consisting of three components: $\varOmega _{C}$ (khaki) represents the subdomain consisting of only coarse grains; $\varOmega _{\textit{CF}}$ (brown with khaki spots) represents the subdomain consisting of fine and coarse grains; and $\varOmega _{F}$ (brown) represents the subdomain consisting of only fine grains. Light blues representing porewater fill the spaces between khaki circles representing coarse grains in the coarse domain $\varOmega _C$ . The brown spaces representing fine grains fill the spaces between the khaki circles representing coarse grains in the coarse domain $\varOmega _C'$ . The domain $\varOmega _F$ consists of only fine grains, occupying the spaces between coarse grains in $\varOmega _C'$ , as well as having its own part of the domain in the large unit cell $\varOmega ''$ . Light blues representing pore water fill the spaces between brown circles representing fine grains in $\varOmega _F$ . The wide variance in the size of brown circles corresponds to the wide variance often seen in the size of fine grains. We indicate a length scale for each domain: $l'''$ corresponds to the system domain scale, $l''$ corresponds to the length scale for compositional fluctuations, $l'$ corresponds to the coarse grain size and $l$ represents the fine grain size.

Figure 4.

Diagram depiction of our conceptualised model domain. The domain $\varOmega '''$ with attached thickness $l'''$ is the overall system scale for which we are interested in estimating permeability. The inset subdomain $\varOmega ''$ consists of three component subdomains: $\varOmega _C$ (khaki) with entirely coarse grains, $\varOmega _{\textit{CF}}$ (brown with khaki hatches) with fine grains in the pores of coarse grains and $\varOmega _F$ (brown) with entirely fine grains. Variations in these subdomains happen at length scale $l''$ . Both subdomains $\varOmega _C$ and $\varOmega _C'$ insets have khaki coarse grains with their sizes $l'$ defining the length scale, but $\varOmega _C$ is filled with porewater (light blue) and $\varOmega _C'$ has fines (brown) within its pores. The subdomain $\varOmega _F$ has brown fine grains defining its length scale $l$ and the pores of $\varOmega _F$ are filled with porewater (light blue). The subdomain $\varOmega _F$ has the smallest length scale, indicated with the highest associated powers of $\epsilon$ .

The fine-component fraction, $c$ , is the volume fraction of the grain-size distribution belonging to the fine mode. For our idealised experiments, we determine $c$ as we create the sample. For our bimodal tills, we identify the fine mode from the coarse mode by identifying the grain size $d_{s}$ corresponding to the lowest local minimum in the grain-size distribution. The fine-component fraction, $c$ , corresponds to the total area of the grain-size distribution below $d_{s}$ . Then, $\phi _T$ is the bulk porosity of the granular medium and is measured directly in the laboratory. The bulk porosity, $\phi _T$ , represents the total volume through which water flows per unit volume of bulk material.

In our notation, $\phi _i(\varOmega _{\kern-1pt j})$ represents the porosity of the component $i$ in the $\varOmega _{\kern-1pt j}$ subdomain. We take $\phi _{C}(\varOmega _{C})$ to be the porosity of the coarse grains in $\varOmega _{C}$ and $\phi _{F}(\varOmega _{F})$ to be the porosity of the fine grains in $\varOmega _{F}$ . Since $\varOmega _{C}$ only has coarse grains and $\varOmega _{F}$ only has fine grains, $\phi _{F}(\varOmega _{C})$ and $\phi _{C}(\varOmega _{F})$ are not needed. However, $\varOmega _{\textit{CF}}$ has both fine and coarse grains that contribute to its bulk porosity, and these porosities may be distinct from their counterparts in other subdomains. We take $\phi _{C}(\varOmega _{\textit{CF}})$ to be the porosity of the coarse grains in $\varOmega _{\textit{CF}}$ . We take $\phi _{F}(\varOmega _{\textit{CF}})$ to be the porosity of the fine grains that fill the pores between the coarse grains in $\varOmega _{\textit{CF}}$ .

Within a representative volume element, the portion of the domain $\varOmega ''$ composed of $\varOmega _{C}$ , $\varOmega _{F}$ and $\varOmega _{\textit{CF}}$ varies. We take these subdomain fractions to be $\varPhi _{C}$ , $\varPhi _{F}$ and $\varPhi _{\textit{CF}}$ , respectively, so that $\varPhi _{C} + \varPhi _{F} + \varPhi _{\textit{CF}} = 1$ .

The total volume of fine grains $v_{F,T}$ and the total volume of coarse grains $v_{C,T}$ in a unit volume provides two additional constraints. We use the corresponding notation convention for $v$ as we used for $\phi$ , so that $v_i(\varOmega _{\kern-1pt j})$ represents the volume of the component $i$ in the $\varOmega _{\kern-1pt j}$ subdomain. Contributions to the total volume of coarse grains in a representative volume element come from the $\varOmega _C$ subdomain or from the $\varOmega _{\textit{CF}}$ subdomain. Contributions to the total volume of fine grains in a representative volume element come from the subdomain $\varOmega _F$ or from the pores of the coarse grains in the subdomain $\varOmega _{\textit{CF}}$ . Using our notation, this comes out as follows:

(3.1) \begin{align} v_{C,T} &= \varPhi _Cv_C(\varOmega _C) + \varPhi _{\textit{CF}}v_C(\varOmega _{\textit{CF}}), \nonumber \\[6pt]v_{F,T} &= \varPhi _F v_F(\varOmega _F) + \varPhi _{\textit{CF}} \phi _C(\varOmega _{\textit{CF}}) v_F(\varOmega _{\textit{CF}}) \, . \end{align}

From our information about the bulk, $1-\phi _T$ is the total bulk volume fraction of solids, meaning the total volume fraction of fines is $v_{F,T} = c(1-\phi _T)$ and the total volume of coarse grains is $v_{C,T} = (1-c)(1-\phi _T)$ . Because for a representative volume element, all $v_i(\varOmega _{\kern-1pt j}) = 1 - \phi _i(\varOmega _{\kern-1pt j})$ ,

(3.2) \begin{align} c(1-\phi _T) &= \varPhi _{F}(1-\phi _{F}(\varOmega _{F})) + \varPhi _{\textit{CF}}\phi _{C}(\varOmega _{\textit{CF}})(1-\phi _{F}(\varOmega _{\textit{CF}})), \nonumber \\[6pt](1-c)(1-\phi _T) &= \varPhi _{\textit{CF}}(1-\phi _{C}(\varOmega _{\textit{CF}})) + \varPhi _{C}(1-\phi _{C}(\varOmega _{C})) \, . \end{align}

We assume that particle sorting does not occur so that a representative volume element can have either a $\varOmega _{F}$ subdomain or a $\varOmega _{C}$ subdomain, but not both. Mathematically, this corresponds to the constraint $\varPhi _{F}\varPhi _{C}=0$ . The assumption is motivated by observations of ring shear experiments (e.g. Hooyer & Iverson Reference Hooyer and Iverson2000) under high normal stress in which shear strain results in mixing of compositions rather than separation.

Our system for six conceptual properties ( $\varPhi _{C}, \varPhi _{F}, \phi _{F}(\varOmega _{\textit{CF}}), \phi _{C}(\varOmega _{\textit{CF}}), \phi _{F}(\varOmega _{F}), \phi _{C}(\varOmega _{C})$ ) still remains underconstrained. For our study, we assume $\phi _{F}(\varOmega _{F})$ and $\phi _{C}(\varOmega _{C})$ to be material constants, and we assume that $\varPhi _C=\varPhi _C(c)$ and $\varPhi _F=\varPhi _F(c)$ are linear in the fine-component fraction. During calibration, we lower bound all porosities using the densest close packings of spheres, $\phi = 0.26$ . We upper bound all porosities using the gravitational packing corresponding to the material. Once the functions $\varPhi _C(c), \varPhi _F(c)$ are calibrated, our system for conceptual properties is fully constrained.

3.2. Sequential homogenisation of Stokes and seepage flow

We solve the momentum equation in the limit of negligible inertia assuming that the pore fluid is incompressible and satisfies no-slip on solid surfaces. This limit corresponds to the following incompressible Stokes problem in index notation:

(3.3) \begin{align} \begin{aligned} \partial _i u_i &= 0, \\[4pt]\partial _i p &=\mu \partial _{\kern-1pt j} \partial _{\kern-1pt j} u_i, \\[4pt]u_i &= 0 \quad \textrm{on} \, \varGamma \, , \end{aligned} \end{align}

where $\varGamma$ is the boundary of any solid surfaces in the computational domain.

From here, we follow the sequential homogenisation procedure reflective of Mei & Auriault (Reference Mei and Auriault1989), establishing two distinct types of coarse-graining. At the smallest scale where the fluid–solid pore geometry is explicitly resolved, we homogenise the Stokes problem to obtain an effective Darcy permeability. At all larger scales, the pore space has already been represented by an effective permeability, so we perform consecutive homogenisation steps of seepage problems through a heterogeneous Darcy medium to obtain bulk permeabilities applicable at larger scales.

Thus, we solve Stokes cell problems that do not always begin at the same physical scale, beginning when the fine-grain scale is resolved as in $\varOmega _F$ , or beginning at coarser grain-scales as in $\varOmega _C$ . Once the first Stokes homogenisation has been performed, all subsequent steps use the same seepage homogenisation procedure iteratively. Rather than carrying out a single perturbation expansion over all four length scales of the system problem, we illustrate the two homogenisation steps as a sequence of adjacent-scale problems. At each step, $\epsilon$ denotes the ratio between two neighbouring length scales.

3.2.1. Homogenisation of the Stokes problem

Let $\varOmega$ be the periodic representative cell at the smallest scale for which the fluid–solid pore geometry is resolved and let $\varGamma$ be the solid boundaries of grains. We define $\epsilon$ such that $l \sim \epsilon l'$ . For coarse-graining through this $l$ -scale, we first assume that the characteristic pressure gradient operates at the adjacent scale $l'$ , suggesting the following non-dimensionalisation:

(3.4) \begin{align} x_i \sim lx_i^* , \quad p \sim Pp^* ,\quad u_i \sim \frac {Pl^2}{\mu l'}u_i^* . \end{align}

After removing the $^*$ superscript, our non-dimensional equations are:

(3.5) \begin{align} \begin{aligned} \partial _i u_i &= 0, \\ \partial _i p &= \frac {l}{l'}\partial _i \partial _i u_{\kern-1pt j} = \epsilon \partial _{\kern-1pt j} \partial _{\kern-1pt j} u_i, \\ u_i &= 0 \quad \textrm{on} \, \varGamma . \end{aligned} \end{align}

We define a corresponding spatial coordinate such that spatial variations in $\boldsymbol{x} := x_i$ are at the finest-scale and $\boldsymbol{x'} := x'_i$ are at the adjacent, coarse-grained scale:

(3.6) \begin{align} \boldsymbol{x'} := x'_i \sim \frac {l}{l'} x_i =\epsilon x_i . \end{align}

This means that the gradients at the different scales are

(3.7) \begin{align} \partial _i' := \frac {\partial }{\partial x'_i} \sim \frac {\partial }{\partial (\epsilon x_i)} = \frac {1}{\epsilon } \partial _i \, . \end{align}

We then conduct a perturbation expansion in the pressure and velocity fields to match the Stokes momentum equation to leading order:

(3.8) \begin{align} p &= p^{[0]}+\epsilon p^{[1]}+\textrm{h.o.t.} , \\[-12pt] \nonumber \end{align}

(3.9) \begin{align} u_i &= u_i^{[0]}+\epsilon u_i^{[1]}+\textrm{h.o.t.} . \\[10pt] \nonumber \end{align}

For our perturbed Stokes problem, the incompressibility condition $\partial _i u_i = 0$ expands to

(3.10) \begin{align} ( \partial _i + \epsilon \partial '_i ) & \left(u_i^{[0]} + \epsilon u_i^{[1]} \right) = 0, \\[-12pt] \nonumber \end{align}

(3.11) \begin{align} O( \epsilon ^{0} ) &: \partial _i u_i^{[0]} = 0, \\[-12pt] \nonumber \end{align}

(3.12) \begin{align} O( \epsilon ^{1}) &: \partial _i u_i^{[1]} + \partial '_i u_i^{[0]} = 0 . \\[10pt] \nonumber \end{align}

The momentum equation $\partial _i p = \epsilon ^3 \partial _i \partial _i u_{\kern-1pt j}$ expands to

(3.13) \begin{align} \left(\partial _i + \epsilon \partial '_i \right) \big(p^{[0]} + \epsilon p^{[1]} \big) = \epsilon \left( \partial _i + \epsilon \partial '_i \right)^2 \left(u_{\kern-1pt j}^{[0]} + \epsilon u_{\kern-1pt j}^{[1]} \right) , \end{align}

with

(3.14) \begin{align} O(\epsilon ^{0}) &: \partial _i p^{[0]} = 0, \\[-12pt] \nonumber \end{align}

(3.15) \begin{align} O(\epsilon ^{1}) &: \partial '_i p^{[0]} + \partial _i p^{[1]} = \partial _{\kern-1pt j} \partial _{\kern-1pt j} u_i^{[0]}. \\[10pt] \nonumber \end{align}

The no-slip boundary condition at the grain boundaries applies to all perturbation terms of the velocity $u_i^{[0]}\rvert _{\varGamma }= u_i^{[1]}|_{\varGamma } = 0$ . However, the assumption of periodicity on $\varOmega$ imposes a periodic boundary condition on $u_i^{[0]}, u_i^{[1]}$ and $p^{[0]}, p^{[1]}$ on cell boundaries $\partial \varOmega \setminus \varGamma$ .

Our (3.14) implies that $p^{[0]}$ only varies in $\boldsymbol{x'}$ , i.e. $p^{[0]}(\boldsymbol{x'})$ . If we view $\partial _i$ and $\partial _l\partial _l$ in (3.15) as linear operators acting on $p^{[1]}$ and $u_i^{[0]}$ , it follows that there should exist some solution pair ( $p^{[1]}$ , $u_i^{[0]}$ ) that can be expressed directly as a linear transformation of $\partial _i' p^{[0]}$ , i.e.

(3.16) \begin{align} p^{[1]} &= S_i\partial '_i p^{[0]} , \\[-12pt] \nonumber\end{align}

(3.17) \begin{align} u_i^{[0]} &= -k_{\textit{ij}}\partial '_{\kern-1pt j} p^{[0]} , \\[10pt] \nonumber \end{align}

where tensor $k_{\textit{ij}}(\boldsymbol{x}, \boldsymbol{x'})$ and vector $S_i(\boldsymbol{x}, \boldsymbol{x'})$ are functions of all spatial scales. Inserting (3.16) and (3.17) back into the momentum balance, (3.15), and the incompressibility condition, (3.11), we can establish a set of equations for $k_{\textit{ij}}$ and $S_i$ :

(3.18) \begin{align} \begin{aligned} \partial _i k_{\textit{ij}} &= 0, \\[3pt]\delta _{\textit{ij}} + \partial _i S_{\kern-1pt j} &= \partial _l \partial _l k_{\textit{ij}} ,\\[3pt]k_{\textit{ij}} = 0 &\quad \textrm{on} \, \varGamma, \\[3pt]k_{\textit{ij}}, S_{\kern-1pt j} \quad &\textrm{periodic over} \, \varOmega. \end{aligned} \end{align}

With this solution, we can now perform an averaging procedure over this smallest scale cell $\varOmega$ . The averaging operator $\left \langle f \right \rangle _{{\varOmega }}$ over $\varOmega$ appears as

(3.19) \begin{align} \left \langle f \right \rangle _{{\varOmega }}(\boldsymbol{x'}) = \frac {1}{|\varOmega |}\int _\varOmega \text{d}V_{\boldsymbol{x}} f(\boldsymbol{x}, \boldsymbol{x'}) \, , \end{align}

where $\left \langle f \right \rangle _{{\varOmega }}$ is only a function of the larger length scales. The upscaling of the velocity $u_i^{[0]}$ only requires the upscaling of the tensor $k_{\textit{ij}}(\boldsymbol{x}, \boldsymbol{x'})$ into a slowly varying tensor $k_{\textit{ij}}' := \left \langle k \right \rangle _{{\varOmega ,\textit{ij}}}( \boldsymbol{x'})$ . This upscaled tensor is the homogenised Darcy permeability for the next consecutive scale and is symmetric positive definite (Mei & Auriault Reference Mei and Auriault1989). Applying the averaging procedure to the incompressibility condition (3.12) demonstrates that the Darcy permeability tensor solves the corresponding Darcy’s law:

(3.20) \begin{align} \frac {1}{|\varOmega |}\int _\varOmega \text{d}V_{\boldsymbol{x}} \left(\partial _i u_i^{[1]} + \partial '_i u_i^{[0]} \right) &= 0 , \\[-28pt] \nonumber\end{align}

(3.21) \begin{align} \frac {1}{|\varOmega |}\int _{\partial \varOmega } \text{d}S_{\boldsymbol{x}}n_i u_i^{[1]} + \partial '_i \left \langle u_i^{[0]} \right \rangle _{{\varOmega }} &= 0, \\[-28pt] \nonumber\end{align}

(3.22) \begin{align} \partial '_i \left \langle u_i^{[0]} \right \rangle _{{\varOmega }} &= 0 , \\[-28pt] \nonumber \end{align}

(3.23) \begin{align} -\partial '_i \left (k'_{\textit{ij}}\partial '_{\kern-1pt j}p^{[0]}\right ) &= 0 \, . \\[0pt] \nonumber\end{align}

3.2.2. Homogenisation of the Poisson-like seepage problem

At larger scales, the pore geometry is no longer resolved explicitly. Instead, the medium is described by a spatially variable Darcy permeability tensor $k_{\textit{ij}}(\boldsymbol{x},\boldsymbol{x'})$ . Let $\varOmega$ now represent the porous domain for the seepage problem and now let $\varGamma$ represent the impermeable internal boundaries of $\varOmega$ . Starting from Darcy’s law that was the result of the homogenisation of the Stokes problem, we can recover the general equation without imposed length scales:

(3.24) \begin{align} \partial _i \left (k_{\textit{ij}} \partial _{\kern-1pt j} p\right ) = 0. \end{align}

The homogenisation of the seepage problem at this scale is formulated as large scale pressure gradients over $l' \sim \epsilon ^{-1}l$ are balanced by small scale variability in mass transport through pores at scale $l$ . We expand the pressure again over one order in $\epsilon$ and use the same convention:

(3.25) \begin{align} p &= p^{[0]}+\epsilon p^{[1]}+\textrm{h.o.t.}\, , \end{align}

with $p^{[0]}=p^{[0]}(\boldsymbol{x'})$ . Substituting the pressure expansion into (3.24), we find that

(3.26) \begin{align} \big(\partial _i + \epsilon \partial '_i \big) \left ( k_{\textit{ij}} \big(\partial _{\kern-1pt j}+\epsilon \partial '_{\kern-1pt j} \big) \big(p^{[0]} + \epsilon p^{[1]} \big) \right )= 0 \, , \end{align}

which, to leading order, has

(3.27) \begin{align} O(\epsilon ^1): \partial _i \left ( k_{\textit{ij}}\big(\partial '_{\kern-1pt j} p^{[0]} + \partial _{\kern-1pt j} p^{[1]} \big)\right )= 0 \, . \end{align}

Once again, the $\partial$ gradient can be treated as a linear operator similar to the homogenisation of the Stokes problem. As a result, (3.27) can be satisfied by a solution where $p^{[1]}$ is related to $\partial '_{\kern-1pt j} p^{[0]}$ by a linear transformation, i.e.

(3.28) \begin{align} p^{[1]} = S_{\kern-1pt j}\left (\partial '_{\kern-1pt j} p^{[0]}\right ) \, . \end{align}

Substituting this expression for $p^{[1]}$ into (3.27),

(3.29) \begin{align} \partial _i \left [k_{\textit{ij}} \left ((\delta _{jl} + \partial _{\kern-1pt j} S_l)\partial '_l p^{[0]}\right )\right ] = 0 \, \quad \forall \partial '_l p^{[0]} \, , \end{align}

we obtain the following vector Poisson equation for the value of $S_{\kern-1pt j}$ over the domain $\varOmega '$ :

(3.30) \begin{align} \begin{aligned} \partial _i (k_{\textit{ij}} \delta _{jl}) &= -\partial _i (k_{\textit{ij}} \partial _{\kern-1pt j} S_l), \\[4pt]n_i k_{\textit{ij}}\left (\delta _{jl}+\partial _{\kern-1pt j} S_l\right )&=0 \quad \textrm{on} \, \varGamma, \\[4pt]S_l \quad &\textrm{periodic over} \, \varOmega , \end{aligned} \end{align}

with a zero normal flux boundary condition corresponding to that used for seepage flow.

With this solution, we can now perform our averaging procedure on the seepage Darcy flux within (3.27) over the seepage $\varOmega$ ,

(3.31) \begin{align} q_i^{[0]} &= \left \langle k_{\textit{ij}}\left(\partial '_{\kern-1pt j} p^{[0]} + \partial _{\kern-1pt j} p^{[1]} \right) \right \rangle _{{\varOmega }} \\[-12pt] \nonumber \end{align}

(3.32) \begin{align} &= \left \langle k_{\textit{ij}}(\delta _{jl} + \partial _{\kern-1pt j} S_l)\left(\partial '_l p^{[0]} \right) \right \rangle _{{\varOmega }} \\[-12pt] \nonumber \end{align}

(3.33) \begin{align} &= \left \langle k_{\textit{ij}} \left(\delta _{jl} + \partial '_{\kern-1pt j} S'_l \right) \right \rangle _{{\varOmega }} \left(\partial '_l p^{[0]} \right)\!, \\[10pt] \nonumber \end{align}

so that the homogenised permeability relevant for the next larger scale is

(3.34) \begin{align} k'_{\textit{ij}} = \left \langle k_{\textit{ij}} \big(\delta _{jl} + \partial _{\kern-1pt j} S_l \big) \right \rangle _{{\varOmega }}\!, \end{align}

and is once again symmetric positive definite (Mei & Auriault Reference Mei and Auriault1989). Averaging (3.26) once again demonstrates that the homogenised permeability tensor solves the corresponding Darcy’s law for the next scale:

(3.35) \begin{align} \partial '_i \left (k_{\textit{ij}} \left(\partial '_{\kern-1pt j} p^{[0]} \right) \right ) = 0\, . \end{align}

Because the output of this homogenisation is also a Darcy permeability tensor, the seepage homogenisation can be applied recursively.

3.2.3. Recursive application to subglacial tills

We now apply these one-step homogenisation procedures to the four length scales of the till problem, from the fine grains between coarse grains to the system scale. In this procedure, we label $\varOmega$ , $\varOmega '$ , $\varOmega ''$ and $\varOmega '''$ as the periodic representative unit cells for each length scale. We define $\epsilon$ such that $l \sim \epsilon l' \sim \epsilon ^2 l'' \sim \epsilon ^3 l'''$ . For subglacial tills, the fine clay scale is $l \sim 10^{-6}$ m and the coarse sand scale is $l' \sim 10^{-3}$ m. For each of our length scales, we define a corresponding spatial coordinate such that spatial variations in $\boldsymbol{x} := x_i$ are at the finest-scale and $\boldsymbol{x'''} := x'''_i$ are at the coarsest-scale:

(3.36) \begin{align} \boldsymbol{x'} := x'_i \sim \frac {l}{l'} x_i =\epsilon x_i , \quad \boldsymbol{x''}:= x''_i= \frac {l}{l''} x_i = \epsilon ^2 x_i , \quad \boldsymbol{x'''}:= x'''_i= \frac {l}{l'''} x_i = \epsilon ^3 x_i\, . \end{align}

These three scales of unit cells apply to all three directions for the purposes of our analysis, meaning that the three vectors corresponding to spatial coordinates can be scaled equally in all directions.

We illustrate the upscaling from the fine grains between coarse grains, through compositional variability and ending at the system scale. If there are regions where the first resolved pore geometry occurs at $\varOmega _C$ rather than at $\varOmega _F$ , then the Stokes homogenisation will start at the coarse-grain scale directly to then be input into the final seepage homogenisation over $\varOmega ''$ . For this example, we go from

(3.37) \begin{align} \varOmega _F \xrightarrow {\langle \rangle } \varOmega _C' \xrightarrow {\langle \rangle } \varOmega '' \xrightarrow {\langle \rangle } \varOmega '''\, . \end{align}

When fine-scale pore geometry is resolved, the first step is the Stokes homogenisation over $\varOmega$ , an operation we refer to as $\mathcal{S}$ . We then undergo consecutive homogenisation over the Poisson-like seepage problem, an operation we refer to as $\mathcal{P}$ . This means that in our illustrative example, the sequential homogenisation is

(3.38) \begin{align} \varOmega _F \xrightarrow {\mathcal{S}} \varOmega _C' \xrightarrow {\mathcal{P}} \varOmega '' \xrightarrow {\mathcal{P}} \varOmega '''\, . \end{align}

We detail a full version of the sequential homogenisation across the length scales in the supplementary material available at https://doi.org/10.1017/jfm.2026.11939, but the results are identical to those from using the recursive procedure outlined previously. The homogenisation procedure leads to a homogenised fine-scale permeability:

(3.39) \begin{align} k'_{\textit{ij}}(\boldsymbol{x'},\boldsymbol{x''},\boldsymbol{x'''})=\left \langle k_{\textit{ij}} \right \rangle _{{\varOmega }} \, \end{align}

from the Stokes homogenisation over $\varOmega _F$ ; a homogenised coarse-scale permeability:

(3.40) \begin{align} k''_{\textit{ij}}(\boldsymbol{x''},\boldsymbol{x'''})=\left \langle \left \langle k_{\textit{ij}} \right \rangle _{{\varOmega }}\left (\delta _{jl}+\partial '_{\kern-1pt j} S'_l\right ) \right \rangle _{{\varOmega ' }} \, \end{align}

from a first seepage homogenisation; and finally, a homogenised compositional-scale permeability:

(3.41) \begin{align} k'''_{\textit{ij}}(\boldsymbol{x'''}) = \left \langle k''_{\textit{ij}}\left (\delta _{jl}+\partial ''_{\kern-1pt j} S''_l\right ) \right \rangle _{{\varOmega '' }} \, \end{align}

from a second seepage homogenisation to be used at the system scale. The final remaining closure from the final seepage homogenisation is the Darcy flow law:

(3.42) \begin{align} \partial '''_i \left (k'''_{\textit{ij}} \partial '''_{\kern-1pt j} p^{[0]} \right ) = 0\, . \end{align}

For comparison with one-dimensional laboratory measurements and scalar relations such as the Kozeny–Carman relation, we report the tensor component in the direction of the imposed pressure gradient, $k = k'''_{ll}$ . When reporting distributions of permeabilities over different geometries, we treat each direction of a geometry realisation independently.

3.3. Numerical solution

We are interested in how different grain-size distributions form different granular packings and thereby different pore geometries to affect the relationship between permeability and porosity. Numerical solution of the homogenised equations from § 3.2 then becomes a computationally efficient approach to accurately see the effects of pore geometry. There are two types of equations that we need to solve: a problem that has the form of a forced Stokes problem at the scale of the smallest pores (3.18) and a Poisson-like seepage problem (3.30) at each consecutive scale that follows.

3.3.1. Fine-scale Stokes problem

We formulate our fine-scale Stokes-like problem using a mixed-element formulation and Nitsche’s method (Juntunen & Stenberg Reference Juntunen and Stenberg2009) to weakly enforce both periodic and Dirichlet boundary conditions. The mixed elements consist of $(k_{\textit{ij}}, S_{l})$ . To simplify notation, we use a shorthand to define operators on the interior and the boundary which splits dot products based on trial function and test function:

(3.43) \begin{align} ( u, v ) = \int _\varOmega \text{d}\boldsymbol{x} uv, \quad [n_i u, v] = \int _\varOmega \text{d}S(\boldsymbol{x}) n_i uv. \end{align}

Note that with this notation, we will explicitly state where the normal vector fits into the dot product operation when taking boundary integrals. Writing the complete problem weakly and using test functions, $\xi _{\textit{ij}}$ and $\psi _l$ , that correspond to $k_{\textit{ij}}$ and $S_l$ , respectively, we obtain

(3.44) \begin{align} - (\partial _i k_{\textit{ij}},\psi _{{\kern-0.5pt}j}) - (\mu \partial _l \partial _l k_{\textit{ij}} - \delta _{\textit{ij}} - \partial _i S_{\kern-1pt j}, \xi _{\textit{ij}}) = 0, \\[-28pt] \nonumber \end{align}

(3.45) \begin{align} - ( \partial _i k_{\textit{ij}},\psi _{{\kern-0.5pt}j}) - [n_l \mu \partial _l k_{\textit{ij}}, \xi _{\textit{ij}}] + (\mu \partial _l k_{\textit{ij}}, \partial _l \xi _{\textit{ij}})\nonumber \\[5pt]+ [n_i S_{\kern-1pt j}, \xi _{\textit{ij}} ] - (S_{\kern-1pt j}, \partial _i \xi _{\textit{ij}}) = - (\delta _{\textit{ij}}, \xi _{\textit{ij}}). \\[10pt] \nonumber \end{align}

We formulate our Stokes problem to use $H_1(\varOmega )$ elements. We now apply Nitsche’s method to symmetrise the problem and enforce our boundary condition. We collapse our Dirichlet and periodic boundary conditions to the condition that $k_{\textit{ij}}= g_{\textit{ij}}, S_l = h_l$ on $\varGamma$ . After reordering the problem, we add a term to symmetrise the problem and two symmetric regularisation terms to enforce the boundary condition. These regularisations are parametrised by scalar coefficients $\eta _1$ and $\eta _2$

(3.46) \begin{align} (\mu \partial _l k_{\textit{ij}}, \partial _l \xi _{\textit{ij}}) - (\partial _i k_{\textit{ij}}, \psi _{{\kern-0.5pt}j}) - (S_{\kern-1pt j}, \partial _i \xi _{\textit{ij}}) + [n_i S_{\kern-1pt j}, \xi _{\textit{ij}}] - [n_l \mu \partial _l k_{\textit{ij}}, \xi _{\textit{ij}}]\nonumber \\[4pt]\quad +[k_{\textit{ij}} - g_{\textit{ij}}, n_i \psi _{\kern-1pt j}] - [k_{\textit{ij}} - g_{\textit{ij}} , n_l \mu \partial _l \xi _{\textit{ij}}] + \eta _1 [k_{\textit{ij}} - g_{\textit{ij}}, \xi _{\textit{ij}}]= - (\delta _{\textit{ij}}, \xi _{\textit{ij}}) \, \end{align}

and

(3.47) \begin{align} &(\mu \partial _l k_{\textit{ij}}, \partial _l \xi _{\textit{ij}}) - (\partial _i k_{\textit{ij}}, \psi _{{\kern-0.5pt}j}) - (S_{\kern-1pt j}, \partial _i \xi _{\textit{ij}}) + [n_i S_{\kern-1pt j}, \xi _{\textit{ij}}] - [n_l \mu \partial _l k_{\textit{ij}}, \xi _{\textit{ij}}] \nonumber \\[4pt]&\quad + [k_{\textit{ij}}, n_i \psi _{\kern-1pt j}] - [k_{\textit{ij}} , n_l \mu \partial _l \xi _{\textit{ij}}] + \eta _1 [k_{\textit{ij}}, \xi _{\textit{ij}}] \nonumber \\[4pt]&\quad= - (\delta _{\textit{ij}}, \xi _{\textit{ij}}) + [g_{\textit{ij}}, n_i \psi _{\kern-1pt j}] - [g_{\textit{ij}} , n_l \partial _l \xi _{\textit{ij}}] + \eta _1 [g_{\textit{ij}}, \xi _{\textit{ij}}] \, . \end{align}

Equation (3.47) is the complete form of the problem that we solve at the microscale. In three dimensions, $k_{\textit{ij}}$ is a nine-component tensor and $S_l$ is a three-component vector, leading to 12 unknowns at each cell. When fully coupled, this problem can become prohibitively expensive, as the size of the matrix corresponding to the problem then scales with the square of the total number of unknowns. However, we can decouple from any relevant contractions along $j$ , instead solving three separate problems of three components, $k$ , and one component, $S$ , using $[a_{\textit{ij}}, n_i b_{\kern-1pt j}] = [n_i a_{\textit{ij}}, b_{\kern-1pt j}]$ , to obtain

(3.48) \begin{align} &(\mu \partial _l k_{\textit{ij}}, \partial _l \xi _{\textit{ij}}) - (\partial _i k_{\textit{ij}}, \psi _{{\kern-0.5pt}j}) - (S_{\kern-1pt j}, \partial _i \xi _{\textit{ij}}) + [S_{\kern-1pt j}, n_i \xi _{\textit{ij}}] - [n_l \mu \partial _l k_{\textit{ij}}, \xi _{\textit{ij}}]\nonumber \\[4pt]&\quad + [n_i k_{\textit{ij}},\psi _{\kern-1pt j}] - [k_{\textit{ij}} , n_l \mu \partial _l \xi _{\textit{ij}}] + \eta _1 [k_{\textit{ij}}, \xi _{\textit{ij}}] \nonumber \\[4pt]&\quad = - (\delta _{\textit{ij}}, \xi _{\textit{ij}}) + [n_i g_{\textit{ij}}, \psi _{\kern-1pt j}] - [g_{\textit{ij}} , n_l \partial _l \xi _{\textit{ij}}] + \eta _1 [g_{\textit{ij}}, \xi _{\textit{ij}}] . \end{align}

As can be seen in (3.48), the problem can be solved by fixing each component $j$ , and contracting along the $i$ and $l$ indices. Furthermore, the bilinear form corresponding to the left-hand side is symmetric under the exchange of the trial and test functions, as desired when using Nitsche’s method. The right-hand side is a constant defined by test functions and boundary conditions.

3.3.2. Intermediate vector Poisson problem

For the Poisson problems that we solve at intermediate scales like $\varOmega _C'$ and $\varOmega ''$ , we once again apply Nitsche’s method (Juntunen & Stenberg Reference Juntunen and Stenberg2009) to weakly enforce periodic and Dirichlet boundary conditions. We are only solving for the vector $S'_i$ , so no mixed-element formulation is necessary. We use $\zeta _i$ as the test function for this formulation and, for brevity, collapse our periodic and Dirichlet boundary conditions as a function $S'_i = S'_{0,i}$ . We form the weak version of the problem to enable the use of $H_1(\varOmega )$ elements and then we apply Nitsche’s method with the regularisation parameter $\eta _2$ :

(3.49) \begin{align} \big(\partial '_i k'_{\textit{ij}} \big(\delta _{jl} + \partial '_{\kern-1pt j} S'_l \big), \zeta _l \big) &= 0, \\[-10pt] \nonumber \end{align}

(3.50) \begin{align} \big[n_i k'_{\textit{ij}} \partial '_{\kern-1pt j} S'_l, \zeta _l \big] - \big(k'_{\textit{ij}} \partial '_{\kern-1pt j} S'_l, \partial '_i \zeta _l \big)&= - \big(\partial '_i k'_{\textit{ij}} \delta _{jl}, \zeta _l \big) \\[-10pt] \nonumber \end{align}

(3.51) \begin{align} - \big(k'_{\textit{ij}} \partial '_{\kern-1pt j} S'_l, \partial '_i \zeta _l \big) + \big[n_i k'_{\textit{ij}} \partial '_{\kern-1pt j} S'_l, \zeta _l \big] + & \big[S'_l - S'_{0,l}, n_i k'_{\textit{ij}} \partial '_{\kern-1pt j} \zeta _l \big] + \eta _2 \big[S'_l - S'_{0,l}, \zeta _l \big] \\[-10pt] \nonumber \end{align}

(3.52) \begin{align} &= -\big(\partial '_i k'_{\textit{ij}} \delta _{jl}, \zeta _l \big) \\[-10pt] \nonumber \end{align}

(3.53) \begin{align} - \big(k'_{\textit{ij}} \partial '_{\kern-1pt j} S'_l, \partial '_i \zeta _l \big) + \big[n_i k'_{\textit{ij}} \partial '_{\kern-1pt j} S'_l, \zeta _l \big] + & \big[S'_l, n_i k'_{\textit{ij}} \partial '_{\kern-1pt j} \zeta _l \big] + \eta _2 \big[S'_l, \zeta _l \big] \\[-10pt] \nonumber \end{align}

(3.54) \begin{align} &= -\big(\partial '_i k'_{\textit{ij}} \delta _{jl}, \zeta _l \big) + \big[S'_{0,l}, n_i k'_{\textit{ij}} \partial '_{\kern-1pt j} \zeta _l \big] + \eta _2 \big[S'_{0,l}, \zeta _l \big] . \\[10pt] \nonumber \end{align}

3.4. Generation of equivalent pore geometries

A grain-size distribution and porosity do not suffice to determine the pore geometry; the variability in the realised granular packing creates stochasticity in the pore geometry and thereby the permeability. To realise a few possible granular packings for a given grain-size distribution and porosity, we use DEM to simulate their dynamics using a modified version of a large-scale atomic/molecular massively parallel simulator (LAMMPS), known as the LAMMPS improved for general granular and granular heat transfer (LIGGGHTS) (Kloss et al. Reference Kloss, Goniva, Hager, Amberger and Pirker2012). We then use Gmsh (Geuzaine & Remacle Reference Geuzaine and Remacle2009) to generate three-dimensional (3-D) meshes to represent the pores between grains. Finally, we use the FEniCS 3-D finite element method (FEM) framework with PETSc Krylov solvers (Baratta et al. Reference Baratta, Dean, Dokken, Habera, Hale, Richardson, Rognes, Scroggs, Sime and Wells2023) to numerically solve our homogenisation equations for the given pore geometry. We use Taylor-Hood P2P1 elements throughout our numerical procedures because they are known to be stable and consistent for Stokes problems (Bercovier & Pironneau Reference Bercovier and Pironneau1979). We scale our regularisation parameters $\eta _1$ and $\eta _2$ with the inverse of the average cell diameter, as suggested by Juntunen & Stenberg (Reference Juntunen and Stenberg2009).

Homogenisation assumes a unit-cell domain that is repeated periodically (e.g. Bensoussan et al. Reference Bensoussan, Lions and Papanicolaou2011). However, our generated pore geometries do not follow a perfectly random distribution for localisded permeability and are therefore susceptible to correlation effects (Dagan Reference Dagan1989). We limit the effects of correlation by using a large unit cell of sides at least 3 times the size of the median grain diameter.

Using a periodic boundary condition in our DEMs does not create a sufficiently homogeneous pore geometry unaffected by the boundary (e.g. Cheung & O’Sullivan Reference Cheung and O’Sullivan2008). To mitigate the inherent heterogeneity induced by periodic boundaries from the finite extent of the DEM domain, we use the following procedure.

1. (i) Initialise a large domain with thousands of particles sampled from the given grain-size distribution.

2. (ii) Match porosity using either gravitational settling or uniformly growing particle diameters (Kloss et al. Reference Kloss, Goniva, Hager, Amberger and Pirker2012).

3. (iii) Identify and separate a small portion of this domain with hundreds of particles.

4. (iv) Periodically replicate this domain with a small gap with a width of the median grain diameter.

5. (v) Fill in gaps between periodic unit cells with particles sampled from a given grain-size distribution.

We depict an example simulation domain for a periodic unit cell at the fine scale in figure 5.

Figure 5.

Example of a periodic simulation domain for the fine scale as visualised in ParaView (Ahrens et al. Reference Ahrens, Geveci, Law, Hansen and Johnson2005). Inset zooms in on one part of the fine scale to emphasise the smaller spheres of the simulated domain.

3.5. Modelling subdomain variability in the $\varOmega ''$ domains

In our model, the only constraint we have placed on an $\varOmega ''$ representative unit cell from figure 4 is the spatial fractions of each subdomain $\varPhi _C$ , $\varPhi _F$ and $\varPhi _{\textit{CF}}$ for $\varOmega _C$ , $\varOmega _F$ and $\varOmega _{\textit{CF}}$ , respectively. The geometry through which these subdomains fill the unit cell can have important implications for the bulk permeability $k''$ of the $\varOmega ''$ domain. For example, layered media oriented perpendicular to the flow direction will have a bulk permeability very different from that of layered media oriented parallel to the flow direction. Although imagery may be able to help establish the geometry of $\varOmega ''$ empirically, the geometry of $\varOmega ''$ will often be uncertain. There are two options to navigate this uncertainty in our model: (i) modelling $\varOmega ''$ as dominated by a single subdomain; or (ii) treating the geometry of $\varOmega ''$ stochastically to generate a distribution of bulk permeabilities $k''$ .

Modelling $\varOmega ''$ with a single uniform subdomain avoids explicitly representing the geometric variability in $\varOmega ''$ , so this assumption relies substantially on user expertise. For example, assuming $\varOmega ''=\varOmega _{\textit{CF}}$ will rely on user judgments about the amount of scale separation between the coarse grains and fine grains of a bimodal material. Assuming $\varOmega ''=\varOmega _{\textit{CF}}$ for materials rich in either the coarse mode or the fine mode (i.e. $c \sim 0$ or $c \sim 1$ ) can lead to unrealistic porosities $\phi _C(\varOmega _{\textit{CF}})$ or $\phi _F(\varOmega _{\textit{CF}})$ to beget unrealistic estimates of permeability. If $c$ is small, there may be too few fine grains to fill the spaces between coarse grains homogeneously. If $c$ is large, coarse grains may act more like asperities in a compacted medium of fine grains.

Approximating $\varOmega ''$ with a single subdomain type may not suffice for some materials, necessitating a stochastic approach to modelling the geometry of $\varOmega ''$ . We characterise $\varOmega ''$ by a subdomain length scale $l_C'' / l''$ , and by their subdomain volume fraction $\varPhi _C$ and $\varPhi _F$ estimated from their composition (i.e. fine-component fraction $c$ ). With the additional assumption that subdomains are distributed homogeneously throughout the material, we generate a large ensemble of realised $\varOmega ''$ geometries consistent with these $l_C''/l''$ and $\varPhi _C$ descriptors. We emphasise that these synthetic cells do not correspond to distributions of specific, imaged reconstructions in the digital-rock sense. Instead, they are an ensemble of subdomain arrangements conditioned on $\varPhi _C$ and $l_C''/l''$ , intended to show how variability in permeability relates to spatial variability in composition. We solve our intermediate-scale Poisson problem from § 3.3.2 on each realised $\varOmega ''$ geometry to obtain a distribution of bulk permeabilities $k''$ . We use statistics of this $k''$ distribution to estimate both the mean and variance of the bulk permeability at larger scale.

Following our full homogenised model from § 3.1, we only need to solve the intermediate-scale Poisson problem for a two-component medium. In this two-component medium, we use the subscript $C$ to refer to the component with higher permeability $k_C$ and the subscript $F$ to refer to the component with lower permeability $k_F$ . After determining the permeabilities of the sub-compositions $k_C$ and $k_F$ , we use the following procedure leveraging Gaussian fields to generate $\varOmega ''$ geometries.

1. (i) Generate a 3-D Gaussian field at high resolution ( $\gt 10^9$ grid-points). Use a spatial filter of subdomain length scale $l_C''/l''$ to create correlated structures in the field (i.e. smooth the noisy field using a Gaussian filter with kernel size $l_C''/l''$ ).

2. (ii) Compute the $\varPhi _C$ th quantile of this field.

3. (iii) Set all values of the field below the $\varPhi _C$ th quantile to $k_C$ and the remainder of the field to $k_F$ .

We solve our intermediate-scale Poisson problem from § 3.3.2 on the generated $\varOmega ''$ geometry to compute a candidate bulk permeability $k''_i$ . We repeat this procedure to generate a large ensemble of $\varOmega ''$ geometries and thereby compute a distribution of bulk permeabilities $k''$ . We use the statistics of this distribution to compute a central estimate as well as variability bounds.

3.6. Verification and validation of numerical methods

To verify the accuracy of our numerical solver of fine-scale Stokes flows using Nitsche’s method, we perform a convergence analysis of our solver against three different analytical solutions. We first benchmark with a smooth and analytical two-dimensional (2-D) case $\boldsymbol{T1}$ that we derive, with the velocity–pressure pair:

(3.55) \begin{align} \boldsymbol{u_{T1}} = (u_1, u_2) = \left (-e^{x_1}(x_2\cos {x_2} + \sin {x_2}), x_2\sin {x_2}e^{x_1}\right ) \!, \quad p = 2e^{x_1}\sin {x_2}. \end{align}

For our numerical solution, we set the velocities at all boundaries of a square domain $\boldsymbol{u}|_\varGamma = \boldsymbol{u_{T1}}$ . After computing the velocities and pressure of the interior of the domain under varying grid size $h$ , we calculate the $L_2$ -error between our computation and the analytical solution given by (3.55).

For our second benchmark $\boldsymbol{T2}$ , we verify against a 2-D lid-driven cavity (Ghia et al. Reference Ghia, Ghia and Shin1982; Gaskell et al. Reference Gaskell, Gürcan, Savage and Thompson1998). For this simulation, we set the following boundary conditions:

(3.56) \begin{align} \boldsymbol{u_{T2}}(x=-1,y)=\boldsymbol{u_{T2}}(x=1,y)=\boldsymbol{u_{T2}}(x,y=0)=\boldsymbol{0},\quad \boldsymbol{u_{T2}}(x,y=2) = 2, \end{align}

for the boundaries of a square domain. After computing the velocities and pressure for the interior of the domain under varying grid size $h$ , we calculate the $L_2$ residual between our computation and the analytical eigenfunction expansion given by Gaskell et al. (Reference Gaskell, Gürcan, Savage and Thompson1998).

Figure 6.

Numerical verification of our method for solving Stokes equations against analytical solutions with increasing grid size $h$ . Insets include streamlines and vector fields of the flow solution as visualised in ParaView (Ahrens et al. Reference Ahrens, Geveci, Law, Hansen and Johnson2005). Three problems are considered: (a) our contrived example of (3.55); (b) lid-driven cavity; and (c) flow past a Stokes sphere. Pluses $(+)$ indicate residuals in the pressure solution where relevant and shaded circular points $(\bigcirc )$ indicate residuals in the velocity solution. Dotted lines indicate approximate power law fits to error in grid size.

For our third verification case $\boldsymbol{T3}$ , we compute the 3-D Stokes flow past a solid sphere of radius $R$ (Stokes Reference Stokes1851). For this simulation, we mesh a box 3-D domain with a sphere in the centre and we enforce no-slip $\boldsymbol{u_{T3}}=\boldsymbol{0}$ on the boundary of the sphere. Ideally, we could simply enforce the background flow velocity $\boldsymbol{u_{T3}}|_\varGamma = U_\infty$ at the boundaries of the box that give the simulation finite extent. However, the perturbation to the Stokes velocity field by a sphere decays slowly $(r^{-1})$ , where $r$ is the distance from the sphere centre, meaning that our simulation domain would need to be extremely large for the boundary condition $\boldsymbol{u_{T3}}|_\varGamma = U_\infty$ to correspond to the real solution for Stokes flow past a sphere to desired accuracy. Instead, we set the flow velocity of the analytical solution (Stokes Reference Stokes1851) as the boundary condition at the faces of the box to be, in spherical coordinates,

(3.57) \begin{align} & u_r(r,\theta ) = U_\infty \cos \theta \left (1 + \frac {R^3}{2r^3} - \frac {3R}{2r}\right ), \nonumber \\[4pt] & u_\theta (r,\theta )=-U_\infty \sin \theta \left (1-\frac {R^3}{4r^3}-\frac {3R}{4r}\right ), \quad u_\phi =0 ,\end{align}

where $\theta$ corresponds to the angle with respect to the direction of the background flow field $U_\infty$ . After computing the velocities and pressure for the interior of the domain under varying grid size $h$ , we calculate the $L_2$ residual between our computation and the analytical solution given in (3.57) and the following (3.58):

(3.58) \begin{align} p(r) = -\frac {3RU_\infty r\cos \theta }{2r^3} \, . \end{align}

In figure 6, we show the results of our convergence analysis for the three benchmark cases with increasing grid size $h$ . In figure 6(a), we show that for $\boldsymbol{T1}$ , the pressure solution has second-order convergence in $h$ and the velocity solution has third-order convergence in $h$ . In figure 6(b), we show that for $\boldsymbol{T2}$ , the velocity solution only attains first-order convergence in $h$ , a result of the discontinuity of the flow at the top corners of the domain (Ghia et al. Reference Ghia, Ghia and Shin1982). We also plot the third-order convergence in $h$ achieved for the lid-driven cavity a distance of at least $0.1$ away from the discontinuity at the corners. In figure 6(c), we show that for set-up $\boldsymbol{T3}$ , the pressure solution has second-order convergence in $h$ and the velocity solution has third-order convergence in $h$ , demonstrating that our method is robust to discontinuities in the flow field.

Having verified our Stokes solver, we now test the accuracy of the homogenisation technique to approximate Stokes flows with increasing separation of scales. In this analysis, we study a 2-D Stokes flow through arrays of circles with two scale-separated radii and vary the ratio $\epsilon$ of these radii. For the direct numerical simulation of Stokes flow, we enforce a vertical pressure drop through the boundary conditions so that $p(x,y=0) = 1$ and $p(x,y=1) = 0$ . For each $\epsilon$ , we compute an approximate flow field through the domain using our homogenisation method and compare our approximation against the flow field computed using direct numerical simulation.

In figure 7, we compare the approximate field from homogenisation with the simulated flow field through direct numerical simulation under varying size ratios $\epsilon$ . In figure 7(a), shaded circular points indicate the $L_2$ residual of our approximation with respect to the directly simulated flow field. The best fit represented by the dotted line indicates that the $L_2$ residual decays at a rate of $\epsilon ^2$ . The images in figures 7(b) and 7(c) provide a visual illustration of where there are discrepancies between the approximation of the flow field and the computed flow field. As can be seen in figure 7(b), large errors are present where the periodicity in the field of small circles is broken by larger circles. Larger $\epsilon$ corresponds to a larger disruption in the periodicity. As can be seen in figure 7(c), large errors are still present where the periodicity in the field of small circles is broken, but a smaller $\epsilon$ corresponds to a smaller disruption in periodicity and so other errors in the approximation throughout the domain become comparable.

Figure 7.

Numerical validation of our homogenisation method for approximating flow fields operating at separate scales against direct numerical simulations of the flow field. (a) $L_2$ residual of the approximation for the flow field as compared to the direct numerical simulation of the flow field for varying degrees of scale separation through the size ratio $\epsilon$ . (b) Residual of the approximation against the direct numerical simulation for large $\epsilon$ . (c) Residual of the approximation against the direct numerical simulation for small $\epsilon$ . Reds correspond to large residuals and blues correspond to small residuals for both visualisations. Residual fields are visualised in ParaView (Ahrens et al. Reference Ahrens, Geveci, Law, Hansen and Johnson2005).

Figure 7. Long description

Panel a is a log-log plot of L2 velocity residual against epsilon, with an approximately quadratic decline as scale separation increases. Panels b and c map residual magnitude for relatively large and small epsilon. Red regions of larger mismatch surround obstacles at large epsilon, whereas the small-epsilon field is almost entirely dark blue, indicating much lower residuals.

Beyond verification and validation in error, a homogenisation approach offers improvements to numerical tractability of computation and time over direct numerical simulation. We discuss the reduction to problem size, computational speed and the increasing alignment of sequential homogenisation with direct numerical simulation results in detail in the supplementay material.

Figure 8.

Ensemble simulations to quantify deviations from the Kozeny–Carman relation as the variance in a unimodal grain-size distribution increases. (a) Permeabilities against porosities as computed from direct numerical simulations (purple scatter points) and as computed from the Kozeny–Carman relation (dotted line). Darker purples indicate larger $\sigma / d_{50}$ . (b) Best fit grain size for simulations of different $\sigma / d_{50}$ when assuming a Kozeny–Carman relation. A grey trendline plots the best-fit to our expansion for the specific-surface. (c) $\textrm{SSE}$ around the best-fit Kozeny–Carman relation for simulations of different $\sigma / d_{50}$ .

Figure 8. Long description

Panel a plots permeability against porosity. Light points with narrow size distributions cluster near the Kozeny-Carman line, while darker, wider distributions fall below it. Panel b shows the equivalent grain size decreasing sharply for sigma/d50 above about 0.3 and departing from a specific-surface fit. Panel c shows normalized fit error increasing, especially near sigma/d50 of 0.4.