1 Introduction

The packing of mixtures of grains is important in many areas of science and engineering, from materials science (e.g. Yu et al. 1993; Norouzi et al. 2012) and process engineering (Ye and Lin 2006) to pharmaceuticals (Masteau and Thomas 1999; Van Veen et al. 2002). In the earth sciences, applications occur in many areas, including limnology, soil science (Zhang et al. 2011), hydrology (Zhang et al. 2009) and petrophysics (Sakaki and Smits 2015; Nooraiepour et al. 2019). Packing is a major influence on the physical properties of all mixtures of grains of different sizes or shapes. Those properties include porosity (Dias et al. 2004), which is a measure of storage capacity, and permeability (Bernabé and Maineult 2015), which describes fluid flow, but also encompass the ability of a material to transmit and absorb acoustic energy (Leurer and Brown 2008), electric current (Glover 2015) and heat (Wallen et al. 2016).

In our case, we are interested in how granular mixtures represented by clastic rocks have their porosity and permeability controlled by grain size, relative grain size, grain shape (Yu et al. 1993), grain orientation and grain packing (Zhang et al. 1996). We recognise rocks whose macroscopic properties exhibit heterogeneous and anisotropic physical properties as a result of changes in grain packing (e.g. Zhang et al. 1996), and that the packing results from the instantaneous availability of sediment of a particular size or shape during the position, and the dynamics of the depositional process itself (Li et al. 2019).

In real rocks, grain size follows a unimodal or multimodal distribution with different amounts of grains of many different sizes. The mathematical description of such three-dimensional packing is extremely complex. Consequently, we will restrict ourselves to mixtures of two sizes of grain (binary mixtures). Considering the wide application and importance even on mixtures of only two sizes of grain, there has been surprisingly little research carried out on them. Examples of binary mixing porosity models that have been developed include the ideal-packing porosity model of Kamann et al. (2007), the fractional-packing porosity model of Koltermann and Gorelick (1995), the fine packing porosity model of Dias et al. (2004) and the interstitiation/non-cutting replacement porosity (INCR) model of Diyokeugwu and Glover (2019). It should be noted that the model of Diyokeugwu and Glover (2019) contains significant errors in the model equations (predominantly their Eqs. 4a and b), which are corrected in this work.

Permeability can either be calculated directly from the properties of the components of the binary mixture, using some distributive models, or using a model which uses the porosity which arises from the porosity mixing models mentioned above. In the case of distributive models, parallel, perpendicular and random patch distributions are adequately described by weighted arithmetic, harmonic and geometric means (Glover et al. 2006), while some have used variants of the Hashin–Shtrikman bounds or Maxwell’s equation (Maxwell 1881), both of which rely on an analogy between fluid flow and electrical transport in porous media. The distributive model approach suffers from the restriction that the arrangement of the zones consisting of each size of grains needs to be defined, and so a true random mixture cannot be fully modelled.

For standard permeability models which use the porosity of the calculated binary mixture, there are many models which could be cited (see Rashid et al. (2015a, b) for a review), but the most common is that of Kozeny (1927) and Carman (1937). This is often used in the hydrocarbon industry, but, though simple to apply, has been discredited because it provides permeabilities often an order of magnitude too high because the model takes no account of dead-end porosity (Walker and Glover 2010). More recently, various variants of the Revil, Glover, Pezard and Zamora (RGPZ) permeability model (Glover et al. 2006; Rashid et al. 2015a, b) have found increasing success for both conventional and tight, clastic and carbonate rocks.

Figure 1 shows a schematic diagram of a 2D binary mixture of circular discs where the volume fraction of smaller discs increases from zero to unity from the left to the right of the diagram. It is analogous to a similar diagram which may be found in Kamann et al. (2007), but differs because we have included curves for porosity and permeability, which are in accord with the model of porosity and permeability developed later in this work. The figure shows two regimes which represent where each of two processes is dominant. At low volume small grain fractions, the interstitiation process (IP) occurs, where smaller grains fill the interstices between larger grains. At high volume fractions of small grains, the replacement process (RP) occurs, where larger grains replace groups of smaller grains and the pore space between them.

Fig. 1
figure 1

Schematic diagram of the physical changes occurring during the mixing of a binary mixture, together with expected changes in porosity and permeability. The subscripts c and f refer to either porosities or permeabilities arising from the packing of either coarse grains only, or fine grains only, respectively

It would be expected that the porosity arising from binary mixing would depend upon relative grain sizes, grain shapes, packing arrangements and how compact, or ideal, the packing is (Woronow 1986). In earth materials, each of these may vary.

The grain size effect is well recognised as a major control on reservoir quality in clastic reservoirs (e.g. Morad et al. 2010) and depends on the availability of grains of a given size (i.e. the grain size distribution), but is also true if only two grain sizes are present. If at any point during the position both small and large grains are available in unrestricted amount, a complete packing might be expected to occur. In this case, we define complete packing as the full occupancy of all void spaces large enough to take at least one of the smaller particles in the mixture. However, there will be times during deposition where the availability of either the smaller grain or the larger grain is restricted. In the case of the unavailability of larger grains, the void spaces will be filled solely with the smaller grains, and so a single-grained sediment would result locally. In the case of the local unavailability of smaller grains, one would expect for spaces to develop within the sediment that are large enough to contain a small grain but do not do so. We define this situation as partial packing, which leads to a higher porosity than would be expected. It should be noted that complete packing is never guaranteed in the natural world, nor in laboratory experiments. Indeed, partial packing is likely to be more common. In our experiments, we have attempted to approach complete packing as closely as possible. However, it should be remembered that the porosity measurements that have been made in this work represent an upper bound to complete packing porosities. By contrast, however, the theoretical equations developed in this work are for complete packing. Consequently, we might expect our experimental measurements to systematically overestimate porosity rather than underestimate it despite all of the considerations that have been put in place to make the measurements as accurate as possible.

Some shapes of grain, such as a cubic grain, pack extremely efficiently if all grains are oriented in the same direction, but extremely poorly, if alternate grains are rotated through 45°. While real rocks do not exhibit the extreme effects seen in cubic grains, the roughness of real grains affects porosities of samples composed of purely one grain size as well as mixtures.

It should also be noted that there may be other reasons for the end-member porosities to vary. For example, in water-saturated sediments, capillary forces stabilise fine grain packs at higher porosities than coarse packs because the specific surface area of fine grain packs is larger than that of coarser packs.

All of the aforementioned effects control the primary depositional porosity. However, post-depositional or diagenetic processes can modify the primary porosity. These effects are either (1) geomechanical (e.g. Obradors-Prats et al. 2019), such as compaction or the development of fractures at all scales from the microscopic (< 10−5 m) to the megascopic (> 103 m), (2) geochemical (e.g. Wu et al. 2019), such as dissolution and precipitation, or (iii) rarely biological in origin (Li and Li 2014). All of these diagenetic effects can occur globally or locally and can lead to heterogeneity and/or anisotropy in the rock. These diagenetic effects are not considered in this work, but process mixing modelling has been carried out recently (Li et al. 2019).

In summary, this paper contains a theoretical study of binary grain packing supported by significant experimental determinations. The applications of the study are very wide, as indicated in the first paragraph of the introduction. A binary mixture is a significant simplification of many of these applications. Only future studies will be able to determine (1) the extent to which this model is efficacious as it stands in each of the applications, (2) whether modifications need to be made to account for a greater distribution of grain sizes or (3) whether the equations can be altered to take account of secondary processes such as compaction, deformation and cementation. This paper only considers the geometric aspects of how packing affects porosity and permeability, and takes no account of the materials from which the grains are composed.

2 Theory: The Two-Dimensional Case

Let us consider a two-dimensional box of arbitrary shape and dimensions with an area Abox = L2. We will place grains of two sizes into this box and examine their packing characteristics and porosity. Initially, the grains are assumed to be circular, but later we will see that the equations developed in this section do not depend on the shape of the grains. The grains are of two physical sizes of radius ȓ and Ȓ for the small- and large-sized discs, respectively. We define and use two-dimensionless measures of grain size, r and R, which represent the radius of each of the grain sizes scaled by the box size. Consequently, the model presented here is independent of box and grain size and is not subject to edge effects. The scaling relationships are:

$$ \begin{aligned} \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{r} & = rL, \\ \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{R} & = RL. \\ \end{aligned} $$
(1)

This is an important first step as it allows the subsequent equations to be independent of the size of the containing box, as well as not relying on the volume of premixed grain fractions as in Kamann et al. (2007) and Koltermann and Gorelick (1995), and in turn allows the mixed porosity to be presented as a function of parameters describing its final state rather than the state of supposed premixed components.

We will start by first considering the process whereby small discs fit between large discs, leading to a reduction in porosity, which we call the interstitiation process (IP). Subsequently, we will consider the porosity reduction process whereby large discs replace small discs and their associated porosity, which we call the replacement process (RP). There are two types of the replacement process, one which implies the partial cutting of small discs by the large discs and one where the small discs pack onto the surface of the larger discs.

2.1 The Interstitiation Process

Consider the box initially packed with the larger discs with some form of packing as shown in Fig. 2a. The packing may be regular, such as square or hexagonal, but could equally well be random, and even partial. The form of the packing is not defined in any way at this stage other than it results from the addition of N large discs to the box resulting in a structure with a porosity

$$ \phi_{{\rm {c}}} = 1 - N\pi R^{2} , $$
(2)

where ϕc represents a porosity defined by the structure of larger (coarser) discs. Note that this equation does not imply any particular packing, either regular or random, nor does it assume that the box is filled with the larger discs. It is possible for there to be large spaces between other large discs that a disc or discs could occupy, but happen not to be occupied. In other words, the packing need not be complete, but may be partial.

Fig. 2
figure 2

a, b The interstitiation process; a the large disc matrix before filling with smaller discs, b partially completed filling with small discs. c, d The replacement process; c the small disc matrix before replacement, d partially completed replacement with large discs, e close-up of the cutting replacement process showing how some grains are unrealistically cut by the process, f close-up of the non-cutting replacement process

Hence, the maximum number of these large discs which can be fitted into the box at a given porosity ϕc is:

$$ N_{\hbox{max} } = \frac{{L^{2} {\kern 1pt} \left( {1 - \phi_{\text{c}} } \right)}}{{\pi \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{R}^{2} }} = \frac{{{\kern 1pt} \left( {1 - \phi_{\text{c}} } \right)}}{{\pi {\kern 1pt} R^{2} }}. $$
(3)

Here, the value of the porosity arising from the packing of large discs ϕc is completely general. It might arise from some type of regular packing (e.g. cubic or hexagonal), perfect random packing or some inefficient random packing. Consequently, ϕc does not represent a minimum porosity arising from large discs, only the actual porosity arising from a given packing of large discs.

Now let us consider the addition of discs which are small enough to fit between the larger discs (r ≪ R), reducing the porosity. During this process, the number of large grains is held constant and equal to the value given in Eq. (3), N = Nmax. The porosity ϕn resulting from this process of adding n small discs is given by:

$$ \phi_{n} = 1 - N\pi R^{2} - n\pi r^{2} = \phi_{\text{c}} - n\pi {\kern 1pt} r^{2} . $$
(4)

This process is shown in Fig. 2b.

This equation is of limited practical utility because it deals with individual numbers of larger and smaller discs. Consequently, we need to rewrite Eq. (4) as a function of the area fraction of fine discs χAf and of coarse discs χAc. The area fraction of fine discs χAf is given by

$$ \chi_{\text{Af}} = \frac{{A_{\text{f}} }}{{A_{\text{f}} + A_{\text{c}} }} = \frac{{n\pi r^{2} }}{{n\pi {\kern 1pt} r^{2} + N\pi R^{2} }}\quad {\text{and}}\quad \chi_{\text{Ac}} = 1 - \chi_{\text{Af}} , $$
(5)

where Af and Ac are the total areas of fine discs and coarse discs, respectively.

Henceforth, we will express all results as a function of the area fraction of small discs χAf. Rearranging Eq. (5) gives

$$ \frac{n}{N} = \left( {\frac{R}{r}} \right)^{2} \frac{{\chi_{\text{Af}} }}{{\left( {1 - \chi_{\text{Af}} } \right)}}. $$
(6)

There are a maximum number of small discs that can fit between the larger discs. This is given by

$$ n_{\hbox{max} } = \frac{{\left( {1 - N\pi {\kern 1pt} R^{2} } \right)\left( {1 - \phi_{\text{f}} } \right)}}{{\pi r^{2} }}, $$
(7)

where ϕf is the porosity defined by the structure if only the smaller discs were packed together, or when defined at a sufficiently small scale that the definition included no effect of the larger discs. Here, the value of the porosity arising from the packing of small discs ϕf is completely general same way as for ϕc, described earlier. The value of ϕf might arise from some type of regular packing (e.g. cubic or hexagonal), perfect random packing, or some inefficient random packing of just the small discs as they fit between the larger discs. Consequently, ϕf does not represent a minimum porosity arising from small discs, only the actual porosity arising from a given packing of small discs.

Consequently, there is a point at which there is no more space between the larger discs to add smaller discs, and this point defines the minimum porosity that the binary structure can attain. Applying Eq. (4), the minimum porosity is defined by

$$ \phi_{\hbox{min} } = 1 - N_{\hbox{max} } \pi R^{2} - n_{\hbox{max} } \pi {\kern 1pt} r^{2} , $$
(8)

which, upon substituting Eq. (7), gives

$$ \phi_{\hbox{min} } = \phi_{\text{f}} \left( {1 - N_{\hbox{max} } \pi {\kern 1pt} R^{2} } \right) = \phi_{\text{f}} \phi_{\text{c}} , $$
(9)

which is in agreement with the classical result. Consequently, the minimum porosity must lie between zero and unity.

Now the porosity from this process can be expressed as a function of the area fraction of small discs by substituting Eqs. (6) and (3) into Eq. (4) to give

$$ \phi_{n} = \frac{{\phi_{\text{c}} - \chi_{\text{Af}} }}{{1 - \chi_{\text{Af}} }}. $$
(10)

This is an interesting result in that it is independent of individual disc sizes and of their relative sizes. The relationship also contains no information about the shape of the discs. Indeed, all of the control on the porosity exercised by the size, shape and packing of the discs resides solely in the values of ϕc and χAf. The value of ϕc includes the effect of grain size and grain shape associated with the distribution of the larger grains, while χAf includes, via ϕf, the same information for the smaller grains.

2.2 The Replacement Process

Now consider the box initially occupied by n smaller grains (Fig. 2c). The porosity of this finer structure depends upon the number of these smaller grains present. The porosity of this structure is

$$ \phi_{\text{f}} = 1 - n\pi {\kern 1pt} r^{2} . $$
(11)

Once again no particular packing, either regular or random is implied by this equation, nor does it assume that the box is filled with the smaller discs, and once again the packing may be partial, and in general that would be the case. If the packing of these small discs is complete, the maximum number of small discs is

$$ n_{\hbox{max} } = \frac{{{\kern 1pt} \left( {1 - \phi_{\text{f}} } \right)}}{{\pi {\kern 1pt} r^{2} }}. $$
(12)

Now imagine adding large discs to the box. Since the large discs are larger than the small discs, there is generally not sufficient space for the simple addition of the large discs. Instead, space must be ‘cut out’ from the fine grain structure to make the space to introduce each new large grain as shown in Fig. 2d. This process results in the overall reduction of porosity because small grains and their associated porosity are removed from the structure and replaced with solid large grains. The replacement porosity reduction process is not as efficient as interstitiation because in that case only pore space is replaced with solid discs.

The porosity of the resulting structure is

$$ \phi_{N} = 1 - n\pi r^{2} - N\pi {\kern 1pt} R^{2} + N\pi R^{2} \left( {1 - \phi_{\text{f}} } \right) = \phi_{\text{f}} - N\pi R^{2} + N\pi {\kern 1pt} R^{2} \left( {1 - \phi_{\text{f}} } \right), $$
(13)

where the \( + N\pi {\kern 1pt} R^{2} \left( {1 - \phi_{\text{f}} } \right) \) term is an adjustment to account for the fact that the \( - \;N\pi {\kern 1pt} R^{2} \) term accounts twice for the addition of the larger solid material. Resolving Eq. (13) gives

$$ \phi_{N} = \phi_{\rm {f}} \left( {1 - N\pi {\kern 1pt} R^{2} } \right). $$
(14)

Once again, there are a maximum number of large discs that may be added before the larger grains would overlap. This is given by

$$ N_{\hbox{max} } = \frac{{{\kern 1pt} \left( {1 - \phi_{\text{c}} } \right)}}{{\pi R^{2} }}, $$
(15)

which is the same as Eq. (3). Substituting Eq. (15) into Eq. (14) allows the minimum porosity attainable using the replacement process and gives

$$ \phi_{\hbox{min} } = \phi_{\text{f}} \phi_{\text{c}} , $$
(16)

which is the same as attained by the interstitiation process in Eq. (9). This result implies that the end-member structures from both processes have equivalent porosities.

Once again, porosity from the replacement process can be expressed as a function of the area fraction of small discs by substituting Eqs. (6) and (12) into Eq. (14) to give

$$ \phi_{N} = \phi_{\text{f}} \left[ {1 - \left( {1 - \phi_{\text{f}} } \right)\frac{{\left( {1 - \chi_{\rm {Af}} } \right)}}{{\chi_{\rm {Af}} }}} \right]. $$
(17)

However, this equation does not accurately represent the porosity evolving from this process in true binary mixtures. The reason is that the process of ‘cutting-out’ sufficient area of small discs and their associated porosity cuts smaller grains apart. Figure 2e shows an enlargement of Fig. 2d to show the effect. This cutting process leaves the larger grains with solid surface knobbles that increase the matrix fraction and occludes porosity. It is clear, therefore, that Eq. (17) leads to a greater reduction of porosity than would be the case if all discs, large or small always remain intact, which is the realistic case. We call this model the cutting replacement process (CRP).

A modified set of equations is needed to include the added limitation that all discs remain whole. This amounts to recognising all of the knobbles on the large disc peripheries, calculating their areas and removing that area from the porosity and area fraction calculations. This sounds complex, but the calculation is not difficult. First, consider the final binary mixture. It contains N large grains and n small grains. The total areas represented by each of the types of disc are

$$ \begin{aligned} A_{\text{c}} & = N\pi R^{2} , \\ A_{\text{f}} & = \left( {1 - \phi_{\text{f}} } \right)\left( {1 - N\pi {\kern 1pt} R^{2} } \right) = \left( {1 - \phi_{\text{f}} } \right)\phi_{\text{c}} ,\end{aligned} $$
(18)

so the area fraction of small grains is

$$ \chi_{\text{Af}} = \frac{{A_{\text{f}} }}{{A_{\text{f}} + A_{\text{c}} }} = 1 - \frac{{N\pi {\kern 1pt} R^{2} }}{{1 - \phi_{\text{f}} + N\pi R^{2} \phi_{\text{f}} }}, $$
(19)

which can be rearranged to give

$$ N = \frac{{\left( {1 - \phi_{\text{f}} } \right){\kern 1pt} \left( {1 - \chi_{\text{Af}} } \right)}}{{\pi {\kern 1pt} R^{2} {\kern 1pt} \left( {1 - \phi_{\text{f}} + \chi_{\text{Af}} \phi_{\text{f}} } \right)}}, $$
(20)

which, when substituted into Eq. (14) for the replacement process, gives

$$ \phi_{N} = \phi_{\text{f}} \left[ {1 - \left( {1 - \phi_{\text{f}} } \right)\frac{{\left( {1 - \chi_{\text{Af}} } \right)}}{{\left( {1 - \phi_{\text{f}} + \phi_{\text{f}} \chi_{\text{Af}} } \right)}}} \right] = \phi_{\text{f}} \left[ {\frac{{\chi_{\text{Af}} }}{{\left( {1 - \phi_{\text{f}} + \phi_{\text{f}} \chi_{\text{Af}} } \right)}}} \right]. $$
(21)

We call this model the non-cutting replacement process (NRP) which results in a structure where small discs are locally packed around the peripheries of the larger discs, as shown in the enlargement in Fig. 2f. Comparison of Eq. (17) and the first form of Eq. (21) shows the difference between the two approaches. The mixed porosity is in both cases less than the porosity of the pure fine disc phase ϕf, but in the case of the non-cutting replacement process the porosity decreases at a slower rate as the area fraction of small discs decreases when larger discs are added than for the cutting replacement process. This is because the knobbles no longer contribute to the porosity reduction.

The value of χAfcrit = χAf(ϕmin) can be found by setting Eq. (10) equal to Eq. (21) while relabelling χAf(ϕmin) = χAfcrit and then rearranging the result to obtain χAfcrit, which after some manipulation gives

$$ \chi_{\text{Afcrit}} = \frac{{\phi_{\text{c}} - \phi_{\text{c}} \phi_{\text{f}} }}{{1 - \phi_{\text{c}} \phi_{\text{f}} }} = \frac{{\phi_{\text{c}} - \phi_{\hbox{min} } }}{{1 - \phi_{\hbox{min} } }}. $$
(22)

Now substituting Eq. (22) into either Eq. (10) or Eq. (21) gives ϕ = ϕcϕf = ϕmin, confirming that the interstitiation process described by Eq. (10) and the non-cutting replacement process described by Eq. (21) meet at the value (χAfcrit, ϕmin).

In summary, two-dimensional binary mixing can be described by

$$ \begin{aligned} \phi_{n} & = \frac{{\phi_{\text{c}} - \chi_{\text{Af}} }}{{1 - \chi_{\text{Af}} }}\quad {\text{for}}\quad \phi_{n} \ge \phi_{\hbox{min} } \quad {\text{and}}\quad \chi_{\text{Af}} \le \chi_{\text{Afcrit}} , \\ \phi_{N} & = \phi_{\text{f}} \left[ {\frac{{\chi_{\text{Af}} }}{{\left( {1 - \phi_{\text{f}} + \phi_{\text{f}} \chi_{\text{Af}} } \right)}}} \right]\quad {\text{for}}\quad \phi_{n} \ge \phi_{\hbox{min} } \quad {\text{and}}\quad \chi_{\text{Af}} \ge \chi_{\text{Afcrit}} , \\ \phi_{\hbox{min} } & = \phi_{\text{c}} \phi_{\text{f}} = \frac{{\phi_{\text{c}} - \chi_{\text{Afcrit}} {\kern 1pt} }}{{1 - \chi_{\text{Afcrit}} }}\quad {\text{and}}\quad \chi_{\text{Afcrit}} = \frac{{\phi_{\text{c}} - \phi_{\hbox{min} } }}{{1 - \phi_{\hbox{min} } }}. \\ \end{aligned} $$
(23)

These four equations define a triangle of which two of the sides are curved. An example of such curves is represented in Fig. 3 with the interstitiation process curve shown as the continuous red line and the non-cutting replacement process shown by the continuous blue line. The blue dashed line represents the cutting replacement process given by Eq. (17). The grey line represents the least efficient mixing of the two base porosities ϕc and ϕf.

Fig. 3
figure 3

Example implementations of the models presented in Eq. (23) with ϕc = 0.228, representing dense random packing of spheres and ϕf = 0.279 representing uniform hexagonal packing of the smaller discs between the randomly arranged larger discs. The minimum porosity ϕmin  = ϕcϕf  = 0.064 and the critical area fraction χAf = 0.176. IP interstitiation process, CRP cutting replacement process, NRP non-cutting replacement process

It is important to note that these curves are independent of physical or relative disc sizes, or the ratio of the size of the small discs to the large discs. They are also independent of grain shape. Although the equations have been derived using circular discs, they are equally valid for and plane regular or irregular two-dimensional shape of constant area. However, that is not to say that the values of ϕn and ϕN will not be different for, say, circular discs and square plates, it is just that the difference lies in and is controlled by the values of ϕc and ϕf, which are inputs to the equations.

For example, let us suppose that we have a matrix of large discs. These have a maximum random packing porosity of about 0.228 (Hinrichsen et al. 1990), i.e. ϕc = 0.228. Let us take two scenarios. In the first, let us suppose a lower-density random packing of irregularly shaped particles (ϕf = 0.35) fills the porosity left by the larger discs. The porosity of the mixture follows the black (for IP) and the blue curves (for NRP) in Fig. 4a. Suppose now that the fill is replaced by a higher-density regular packing of regularly shaped pentagons, with a consequently lower porosity (ϕf = 0.15). The porosity of the mixture follows the black and red curves in Fig. 4a. Note that the IP curve remains common to both scenarios because the matrix of larger grains has not changed. However, if, for either scenario, the packing of the large grains was decreased such that their porosity was increased to ϕc = 0.35, it would be the grey curve in this figure that would be used with the red and blue curves to define the behaviour of the porosity of the mixture.

Fig. 4
figure 4

a IP and NRP curves for the four scenarios where ϕc  = 0.228, ϕc  = 0.35, ϕf = 0.35 and ϕf  = 0.15, b zero-porosity arrangement of square tiles, c highest porosity (ϕf  = 0.50) arrangement of square tiles, d the solution space shown in pink for porosity where the arrangement of square tiles falls between arrangements (b, c). Note that for χAf  ≤ χAfcrit the solution is unique and follows the black IP curve

Some regular packings have variable porosities that arise from the relative rotations of their particles. This will give rise to a range of behaviours. Let us suppose a dense cubic regular packing of smaller square particles fills the porosity left by a densely packed random arrangement of larger discs (ϕc = 0.228). Now the minimum porosity for a dense cubic regular packing of squares is zero (Fig. 4a), while its maximum porosity is 0.5 (Fig. 4b), depending on how the squares are rotated with respect to each other in the pack, i.e. 0 ≤ ϕf ≤ 0.5. The porosity of the binary mixture can be calculated from Eqs. (22) and results in a range given by the shaded area in Fig. 4d.

3 Theory: The Three-Dimensional Case

All of the steps followed for two-dimensional binary mixing may be calculated for particles in three dimensions. This time considers a three-dimensional box of arbitrary shape and dimensions with a volume Vbox = L3. We will place particles of two sizes into this box and examine their packing characteristics and porosity. Once again the grains are initially assumed to be spherical, but will later be extended to arbitrary shapes. Grains of two physical sizes of radius ȓ and Ȓ for the small- and large-sized particles, respectively, are expressed in a dimensionless form using the scaling relationships in Eq. (1).

Both the interstitiation and both types of replacement process are characterised by carrying out the transformation πr2 → 4πr3/3 and πR2 → 4πR3/3 and defining the volume fraction of small grains as

$$ \chi_{\text{Vf}} = \frac{{V_{f} }}{{V_{f} + V_{c} }} = \frac{{n\pi {\kern 1pt} r^{3} }}{{n\pi {\kern 1pt} r^{3} + N\pi {\kern 1pt} R^{3} }}\quad {\text{and}}\quad \chi_{\text{Vc}} = 1 - \chi_{\text{Vf}} . $$
(24)

This gives an interstitiation porosity

$$ \phi_{n} = 1 - \frac{4}{3}N\pi {\kern 1pt} R^{3} - \frac{4}{3}n\pi {\kern 1pt} r^{3} = \phi_{\text{c}} - \frac{4}{3}n\pi {\kern 1pt} r^{2} , $$
(25)

a replacement porosity

$$ \phi_{N} = 1 - n\pi {\kern 1pt} r^{2} - N\pi R^{2} + N\pi R^{2} \left( {1 - \phi_{\text{f}} } \right) = \phi_{\text{f}} - N\pi R^{2} + N\pi R^{2} \left( {1 - \phi_{\text{f}} } \right), $$
(26)

a volume fraction scaling relationship fraction for small particles χVf of

$$ \frac{n}{N} = \left( {\frac{R}{r}} \right)^{3} \frac{{\chi_{\text{Vf}} }}{{\left( {1 - \chi_{\text{Vf}} } \right)}}, $$
(27)

and a minimum porosity which is the same as for the 2D case:

$$ \phi_{\hbox{min} } = \phi_{\text{f}} \left( {1 - \frac{4}{3}N_{\hbox{max} } \pi {\kern 1pt} R^{3} } \right) = \phi_{\text{f}} \phi_{\text{c}} . $$
(28)

In summary, three-dimensional binary mixing can be described by

$$\begin{aligned} \phi_{\rm{n}} & = \frac{{\phi_{\text{c}} - \chi_{\text{Vf}} }}{{1 - \chi_{\text{Vf}} }}\quad {\text{for}}\quad \phi_{\rm{n}} \ge \phi_{\hbox{min} } \quad {\text{and}}\quad \chi_{\text{Vf}} \le \chi_{\text{Vfcrit}} , \\ \phi_{\rm{N}} & = \phi_{\text{f}} \left[ {\frac{{\chi_{\text{Vf}} }}{{\left( {1 - \phi_{\text{f}} + \phi_{\rm{f}} \chi_{\text{Vf}} } \right)}}} \right]\quad {\text{for}}\quad \phi_{\rm{n}} \ge \phi_{\hbox{min} } \quad {\text{and}}\quad \chi_{\text{Vf}} \ge \chi_{\text{Vfcrit}} , \\ \phi_{\hbox{min} } & = \phi_{\text{c}} \phi_{\text{f}} = \frac{{\phi_{\text{c}} - \chi_{\text{Vfcrit}} {\kern 1pt} }}{{1 - \chi_{\text{Vfcrit}} }}\quad {\text{and}}\quad \chi_{\text{Vfcrit}} = \frac{{\phi_{c} - \phi_{\hbox{min} } }}{{1 - \phi_{\hbox{min} } }}, \\ \end{aligned}$$
(29)

which are the same equations as for the 2D case and contain no explicit information about the size or shape of the particles, nor the efficiency of the packing. Consequently, the operative equations are also invariant under the dimensional transformation between 2D and 3D.

Once again our contention is that information about porosity for either the large particles or the small particles is present in their respective pure phase porosities, and their combination allows the porosity of the binary mixture to be found using Eq. (29). In realistic geological situations, packing is often unlikely to be fully developed. For example, in the process of sedimentation, holes that could accommodate a particle might occur purely because particles of the correct size were not present at the right time. Alternatively, diagenetic process might through dissolution open up pore volumes, while precipitation and compaction processes could reduce them. Each of these primary and secondary porosity processes might lead to values of porosity that differ from the curves defined by Eq. (29), and the extent of the discrepancy could be used to track the progress of such processes.

4 Experimental Verification

An initial set of experiments was carried out as a pilot study and reported by Diyokeugwu and Glover (2019). Here, we report a much more comprehensive suite of experiments that have been carried out on a wider range of binary mixtures and using enhanced measurement and mixing protocols. In each case, the porosity and the permeability of each sample were measured.

Each sample was a mixture of two sizes of glass sphere. The glass spheres were the same type used in Glover and Déry (2010) and are summarised in Table 1. The stock of each size of beads was washed alternately in distilled water, acetone and distilled water to ensure that all traces of glass dust and other impurities had been removed. The stock grains were then oven-dried and stored.

Table 1 Glass beads and sand used to make experimental samples

The apparatus for measuring the porosity and permeability is shown in Fig. 5. Samples were made from two sizes of glass bead in various volume fractions using the dry weight of each fraction to ensure that the grain fraction of the subsequently mixed beads was known accurately, according to the formula

Fig. 5
figure 5

The experimental apparatus developed for the measurement of porosity and permeability in this work

$$ \chi_{\text{Vf}} = \frac{{W_{\text{f}} }}{{W_{\text{f}} + W_{\text{c}} }}. $$
(30)

This approach is accurate, providing that the density of the glass used in the two beads is the same. All of the beads used in this work are from Endecotts Ltd., who have confirmed that they are made from the same grade of soda–lime glass using the same process. Consequently, it is reasonable to assume that Eq. (30) can be used for calculating the experimental volume fraction of fine grains. The glass beads themselves are calibration grade samples, with calibration traceable to the National Physics Laboratory. Nevertheless, each of the samples was analysed by laser diffractometry, the diameter measurements of which are given in Table 1.

The two samples of glass beads were initially mixed in distilled water before being decanted into the sample holder. The fine grid in the sample holder ensures that no beads were washed out of the sample holder. Once in the sample holder the mixture was stirred with a glass rod, the rod was washed with distilled water to ensure that no beads were lost from the mixture. The sample holder was then sealed. The filled sample cell was then purged with air until it was completely dry. The sample holder with the dry binary mixture was then subjected to a regime of agitation, at three different frequencies. Each agitation set consisted of three parts: first, low-frequency manual shaking for 1 min; second, tapping at about 4 Hz using a small rubber hammer driven by a geared 9 V motor for 5 min; and third, applying ultrasound at 45 kHz for 10 min. There were three sets of agitation. The sample cell was adjusted such that the sample occupied the smallest bulk volume possible after each part of each agitation. Finally, the sample was purged with CO2 and then flushed with a 0.1 mol/dm3 solution of deaerated water to saturate it completely. Circulation was carried out until the sample was completely water saturated. The sample was mounted vertically as shown in Fig. 5. A measurement of the electrical conductivity across the sample was carried out in order that the formation factor of the sample could be measured. The agitation process replaced the stirring process that had been used in previous pilot measurements (Diyokeugwu and Glover 2019).

Porosity measurements were calculated from knowledge of the volume of the mixed sample and of the solid beads. The sample volume measurement was made by knowing the cross-sectional area of the sample holder and measuring the length of the core holder occupied by the sample. The volume of the beads was measured using the total weight of the sample, the known density of the beads and assuming that the beads were spherical. A full analysis of the random, systematic and total errors for the porosity measurements is given in “Appendix”.

Permeability measurements were made by timing the passage of a volume of distilled water through the sample using a gravitational head approach. The input fluid was stored in a 5 dm3 translucent Nalgene bottle with an outlet valve. The bottle was placed above the sample holder and calibrated with three horizontal marks 1 cm apart as shown in Fig. 5b. At the start of the experiment, the bottle was arranged such that it was full of 0.1 mol/dm3 NaCl solution to the top mark and that the middle mark was 120 ± 0.1 cm above the centre of the sample. Fluid was flowed until the volume in the input bottle reached the bottom mark. Flow rate was measured gravimetrically by collecting the sample holder effluent in a flask arranged on a laboratory scale. The implementation of a gravimetric method obviates the need to have a range of graduated burettes of different capacities and improves the accuracy of the permeability measurement. The measured mass was logged by computer, which also timed the experiment and made adjustments for the input pressure according to the volume of fluid flowed.

The pressure adjustment requires the previous knowledge of the relationship between volume in the input bottle and the position of the fluid surface with respect to the three reference marks. This calibration was carried out before the experiments.

The density of the process water at the experimental temperature was measured extremely accurately using a graduated flask and the laboratory scale. The length and cross-sectional area of the sample are already known for the porosity measurement. The viscosity of the process water was measured using a high-resolution viscometer. Once again, all of the individual systematic and random errors were quantified and analysed to obtain the random, systematic and total errors for the permeability measurement. This analysis is given in “Appendix”.

Measurements were made flowing fluid through the cell in both directions by reversing the cell physically in the experimental rig. This is extremely easy to arrange as the cell itself is symmetric. At least ten measurements of flow were carried out for each mixture in each flow direction allowing a mean permeability and its associated standard deviation to be calculated. It is this experimental uncertainty that is represented by the error bars for permeability in subsequent figures. It is worthwhile noting that we found no systematic difference on the permeability as a function of flow direction, leading us to assume that the samples we tested were as randomly mixed as possible.

It should be noted that for a few samples there was a trend for the permeability from repeated flow measurements to be systematically progressively lower than those taken earlier in the set of ten measurements. This behaviour occurred typically for samples that had a large contrast in grain size and at small volume fractions of fine particles. We ascribed this behaviour to the movement of the small particles within the mixture in response to the flow of fluid, such that they preferentially blocked flow pathways that were previously open. When this was noticed, the sample was remixed, and the permeability returned to the permeability it had at the beginning of the set of flow measurements. Consequently, for the samples we used the results from 20 flow measurements separated by mixing to calculate the permeability. It should be noted that at no point did we encounter any pressure dependence of permeability, which we took as an indication that the grains were well mixed and the fluids contain no air bubbles.

A total of 165 individual measurements have been made on 12 different binary mixtures of glass beads and sand. A minimum of 12 and a maximum of 17 volume fractions of each mixture were tested. This represents a very wide coverage of the binary space as a function of both volume fraction and size ratio. The resulting porosity and permeability data are given in Table 2.

Table 2 Measured porosity and permeability results

5 Porosity

Figure 6 shows the measured porosities (the data of which are given in Table 2) for binary mixtures of the bead and sand sizes characterised in Table 1 as a function of volume fraction of fine beads or grains χVf. Figure 6 shows that in almost all cases the porosity of the mixtures is approximately in agreement with the theoretical curves. Occasionally, the experimental data points follow the theoretical curves within the marked error bounds; however, there is a tendency for porosity to be overestimated. In other words, where the experimental data points do not fall on the theoretical curves they tend to almost always plot slightly above theoretical curves. We attribute this to not attaining the perfect mixing, despite the complex protocol we use to approach it as closely as possible.

Fig. 6
figure 6figure 6

Measured porosity as a function of volume fraction of fine beads or grains from experimental measurements (symbols with calculated uncertainties). The solid lines represent the results of Eq. (29), which are the theoretical curves for ideal packing developed in this work

Examination of the different parts of Fig. 6 shows that binary mixtures of two sizes of glass beads, where the two diameters differ significantly, tend to produce results which follow the theoretical curves into the critical range of volume fractions the smaller particles, where the interstitiation process and the replacement process change over (e.g. Fig. 6a). By contrast, binary mixtures where the two beads have similar diameters produce experimental porosity values in the critical range that are overestimated (e.g. parts i and j of Fig. 6). We propose that it is difficult to attain perfect mixing for beads of similar diameter, at least with a cell the size of which we use. It might be possible to attain a better mixing in a much larger cell.

Two experiments are included that use sand. The mixture of Ottawa sand and fine sand (e.g. Fig. 6k) shows the same general behaviour as for the glass beads. The experimental data points follow the theoretical curves well, but produce values slightly higher than the theoretical curves. Once again, we attribute this variation to not quite attaining perfect mixing. Figure 6l shows the results from a mixture of the largest glass beads (Sample 12) and Ottawa sand. Once again, the experimental data follow the theoretical curves, but with porosity values slightly higher than the theoretical curves define. In both experiments using natural sand, there was a slightly greater scatter in experimental data. This may be due to the difficulties in obtaining perfect packing when the particles composing the binary mixture are irregular in shape.

Tests were also carried out using a mixture that resulted from mixing in water followed by no shaking. For these samples, it was found that the initial packing of grains provided porosities that were always significantly higher than predicted by our theoretical curves. These mixtures contain voids which are at larger than the smaller of the two grain sizes and form because of a local dearth of the smaller grains during mixing and shaking in our experiments, or during deposition in the case of real rocks. We expect that such pore spaces commonly arise in the deposition of sedimentary rocks and hence expect our equations to represent the lower limit on porosities encountered in nature. Figure 7 shows a plot of how the porosity of one sample evolved after repeated shaking, noting that the final values are those given in Table 2.

Fig. 7
figure 7

The development of the measured porosity for three glass bead mixtures as a function of the amount of shaking they had undergone

6 Permeability

6.1 Introduction

Permeability may be predicted from porosity for granular materials using a very large number of models. Almost all of the models are empirical in nature, having been developed for particular rock types. Some of those conventional models, such as that of Kozeny and Carman (Kozeny 1927; Carman 1937), though not fundamentally empirical are implemented using fitting coefficients and consequently become empirical. The Kozeny–Carman (Kozeny 1927; Carman 1937) model has been used previously by Kamann et al. (2007) in their study of the binary mixing phenomenon. Unfortunately, we cannot reproduce their curves and have consequently concluded that there is an error in Kamann et al.’s implementation of the model leading to values many times greater or less than the Kozeny–Carman relationship produces. Furthermore, it is now known that the Kozeny–Carman relationship is not accurate for most rocks because it does not take into account dead-end and unconnected porosity (Walker and Glover 2010). Consequently, we have taken a different approach to calculating the permeability.

Theoretical permeability curves have been calculated directly using three standard mixing models (arithmetic, harmonic and geometric) as well as with two deterministic models (Kozeny–Carman (KC) and RGPZ Exact (RGPZe)). The standard mixing models mix the known permeabilities of each of the pure phases directly, not using the binary porosity calculated and measured earlier in this work. The deterministic models use the mixed binary porosity calculated using Eq. (29), a constant cementation exponent and an effective grain diameter for the mixture, which has been obtained using four different approaches (cubic, arithmetic, harmonic and geometric). Overall, there were 11 different permeability curves calculated in this paper (three standard mixing models and four for each of the deterministic models) together with, as a comparative reference, the experimental data for the mixtures between bead samples 8 and 12 (porosity data of which are shown in Fig. 6a). Figure 8 contains two parts which show the same data, but on different scales; permeability is shown on a restricted linear scale in Fig. 8a and for the same data shown on a logarithmic scale in Fig. 8b.

Fig. 8
figure 8

Comparison of the 11 models studied in this work for the parameters relevant for mixtures of bead samples 8 and 12 (see Table 1), together with the relevant permeability measurements as reference. a Permeability shown on restricted linear axis, b the same data shown on a log scale. The curves for the three standard mixing models are shown in black. The curves for the RGPZe model with effective grain size calculated with each of the cubic scaling, arithmetic mean, harmonic mean and geometric mean are shown in solid coloured lines, and the curves for the KC model are shown with dashed coloured lines

6.2 Standard Mixing Models

Standard mixing models provide a permeability for a mixture based upon the permeabilities of the end members. One end member is the permeability of a sample composed of only the coarse-grained component, and the other end member is the permeability of the sample composed of only the fine-grained component. These end-member permeabilities have been calculated using the exact RGPZ equation (Glover et al. 2006) in each case. The three standard mixing models in this group are as follows:

  1. 1.

    The weighted arithmetic mean mixing model, which is generally employed to calculate the overall permeability of a set of permeabilities arranged in parallel to the flow direction and for a binary mixture, is given by

    $$k_{\text{arith}} = \chi_{\text{Vf}} k_{\text{f}} + \left( {1 - \chi_{\text{Vf}} } \right)k_{\text{c}} ,$$
    (31)
  2. 2.

    The weighted harmonic mean mixing model, which is generally employed to calculate the overall permeability of a set of permeabilities arranged perpendicular to the flow direction and for a binary mixture, is given by

    $$ k_{\text{harm}} = \left( {\frac{{\chi_{\text{Vf}} }}{{k_{\text{f}} }} + \frac{{\left( {1 - \chi_{\text{Vf}} } \right)}}{{k_{\text{c}} }}} \right)^{ - 1} , $$
    (32)
  3. 3.

    The weighted geometric mean mixing model, which is generally employed to calculate the overall permeability of a set of permeabilities arranged randomly with respect to the flow direction (e.g. Glover et al. 2006) and for a binary mixture, is given by

    $$ k_{\text{geom}} = \exp \left\{ {\chi_{\text{Vf }} \ln \left( {k_{\text{f}} } \right)} \right\} + \left( {1 - \chi_{\text{Vf}} } \right)\ln \left( {k_{\text{c}} } \right). $$
    (33)

Figure 8 shows the standard mixing models as black lines. In each case, the value of permeability intermediate volume fractions fall always between the permeability of the two end members. In addition, difference between the weighted arithmetic mean and weighted harmonic mean mixing models can be considerable when there are approximately equal volume fractions of fine and coarse grains. Clearly, the standard mixing models cannot produce a permeability higher or lower than the range encompassed by the two end-member permeabilities. Immediately, one might question their validity because intuitively one would realise that the permeability of a coarse-grained rock where the pores are filled with very much finer material will be lower than either the finer material or the course material alone. The failure of this type of model is partially due to them not including the effect of a mixed porosity as given by Eq. (29). It is clear that this type of model should never be used to calculate the permeability of binary mixtures. However, it does not preclude the use of the geometric mean approach to model an aggregate of randomly distributed representative elementary volumes of rock each with a different permeability, as considered by Glover et al. (2006).

6.3 Deterministic Models

Two deterministic models have been used in this work. Despite the perceived weaknesses in the Kozeny–Carman, it is widely applied and used in previous work on binary mixing (e.g. Kamann et al. (2007)). Consequently, we have implemented in order that the results from this work can be compared with previous work, using the form

$$ k_{KC} = \frac{{d_{\text{eff}}^{2} \phi_{\text{mix}}^{2} }}{{180\left( {1 - \phi_{\text{mix}} } \right)^{3} }}, $$
(34)

where kKC is the permeability (in D), deff is the effective grain size of the mixture (in m) and ϕmix is the porosity of the mixture.

We have also calculated the permeability using the exact form of the RGPZ equation (Glover et al. 2006; Rashid et al. 2015a, b). The permeability is given by

$$ k_{RGPZe} = \frac{{d_{\text{eff}}^{2} }}{{4am^{2} F\left( {F - 1} \right)^{2} }} = \frac{{d_{\text{eff}}^{2} }}{{4am^{2} \phi_{\text{mix}}^{ - m} \left( {\phi_{\text{mix}}^{ - m} - 1} \right)^{2} }}, $$
(35)

where kRGPZe is the permeability (in m2), deff is the effective grain diameter (in m), ϕ is the porosity (unitless), m is the cementation exponent (unitless), which is theoretically equal to 1.5 for spherical grains (Glover 2009), a is a constant that is thought to be equal to 8/3 for three-dimensional samples composed of quasi-spherical grains and F = ϕm is the formation factor. This equation is analytical, derived from electro-kinetic considerations. Like many equations in electro-kinetics, it assumes that the electrical double layer is thin, which requires the value of the cementation exponent m to be derived from measurements on samples saturated with medium- to high-salinity fluid.

The exact RGPZ equation has three input parameters: the porosity, the cementation exponent and the grain diameter. Its advantage over models that only contain porosity as an input parameter is that the exact RGPZ equation also takes other physical controls on permeability, such as pore connectivity, into account. However, this general advantage is a disadvantage in our case because for a mixture of two grain sizes we do not know which grain size or cementation exponent would be characteristic of any given mixture. The same problem also exists for the Kozeny–Carman equation, which also uses the effective grain size for the mixture.

In our study, we know the porosity, cementation exponent and grain size of each of the pure components, both theoretically and experimentally. In addition, we know the porosity of any mixture from the theoretical curves (Eq. 29) as well as experimental measurements of porosity for 165 different fine particle fractions spread over 12 different binary mixtures. By contrast, we do not know the value of the effective grain diameter to use in Eqs. (3435) for each mixture, and in the case of the RGPZ model, neither do we know the value of the effective cementation exponent of the mixture.

The cementation exponent was measured during the porosity and permeability experimental determinations and found to vary by a very small amount, the average and standard deviation for the pure phases being m = 1.4877 ± 0.0599. The values of cementation exponent vary so little that for the modelling carried out in the rest of this paper, we have assumed that the cementation exponent for all mixes is the simple arithmetic mean of cementation exponents of the end members. Tests carried out using a weighted arithmetic mean were indistinguishable from using the simple arithmetic mean. It is worth noting that while this assumption is that the simple arithmetic mean for calculating cementation exponent of a mixture is valid for all of the mixtures we have tested in this work, there is no guarantee that it would be a valid assumption for all other mixtures, for example a mix of grains and strings.

In all of our experiments, we found that the minimum porosity corresponds very well with the minimum permeability. That arrangement is also borne out by the theoretical development presented in this paper. However, experiments of porosity and permeability in geological materials sometimes show that the minima in porosity and permeability do not accord with each other. The reason for this is that the connectivity of the pore space, which controls the permeability and is expressed through cementation exponent, there is as a function of the volume fraction of small particles in the mixture. In consequence, mixtures with larger porosity may have smaller permeabilities compared to other mixtures. We saw no evidence for this in our work, where our mixtures were isotropic. We hypothesise that it may be more likely to occur when the grains are significantly oblate, and lead to anisotropy of permeability within the mixture.

The mixing of the grain sizes to obtain an effective grain diameter presents a problem which we have attempted to address by using four different mixing models for grain size. These are (1) a cubic model, (2) a volume fraction weighted arithmetic mean, (3) a volume fraction weighted harmonic mean and (4) a volume fraction weighted geometric mean.

The cubic model is simply a cubic scaling of volume fraction between the two elements of the binary mixture to take account of the fact that χVf is a measure of volume fraction in three dimensions, while the grain diameter has one linear dimension. The cubic model is given by

$$ d_{\text{eff}} = d_{\text{c}} - \left( {d_{\text{c}} - d_{\text{f}} } \right)\chi_{\text{Vf}}^{1/3} , $$
(36)

where deff is the calculated grain size, dc and df are the coarse and fine grain diameters, respectively, and χVf is the volume fraction of the fine grains in the mixture.

The volume fraction weighted arithmetic, harmonic and geometric means are given by

$$ d_{\text{eff}} = \chi_{\text{Vf}} d_{\text{f}} + \left( {1 - \chi_{\text{Vf}} } \right)d_{\text{c}} , $$
(37)
$$ d_{\text{eff}} = \left( {\frac{{\chi_{\text{Vf}} }}{{d_{\text{f}} }} + \frac{{\left( {1 - \chi_{\text{Vf}} } \right)}}{{d_{\text{c}} }}} \right)^{ - 1} \quad {\text{and}} $$
(38)
$$ d_{\text{eff}} = \exp \left\{ {\chi_{\text{Vf }} \ln \left( {d_{\text{f}} } \right)} \right\} + \left( {1 - \chi_{\text{Vf}} } \right)\ln \left( {d_{\text{c}} } \right),\quad {\text{respectively}}. $$
(39)

We noted in this work that for the range of χVf (0:1), the results of all of the models given in Eqs. (36) to (39) differed. Figure 9 shows the results of all four models as a function of a linear scaling between extreme points, which for the purposes of this diagram are the diameters of bead samples 8 and 12, respectively (Table 1). It is clear that the arithmetic mean model is simply linear, while the other three models provide effective bead diameters smaller than the arithmetic mean model. The harmonic model provides the largest variation from the arithmetic/linear model, while the geometric mean model occupies the middle ground and the cubic scaling model follows the geometric mean model at small bead diameters and the harmonic mean model at large bead diameters.

Fig. 9
figure 9

Effective bead diameter from each of four different mixing models as a function of the effective bead diameter calculated with a simple linear scaling

Application of the deterministic models does not require the use of end-member permeabilities explicitly, but calculates the permeability directly from Eqs. (34) and (35) using the values of porosity from the mixing models developed in this work (i.e. Eq. (29)), a cementation exponent (either m = 1.5, which is known to be valid for packs of spherical particles (Glover 2015), or m = 1.4877, which we obtained experimentally) and effective grain sizes from each of the models described above.

The deterministic models are represented in Fig. 8 by the solid coloured lines for the RGPZ model and dashed coloured lines for the KC model. It is immediately clear that the minimum porosity present in the porosity models that arises from Eq. (29) is carried over into the permeability model and results in a permeability minimum occurring at the same volume grain fraction. This permeability minimum can, according to the model, provide permeabilities lower than any of the standard mixing models or their end-member permeabilities and consequently conforms better to how we might intuitively imagine the permeability to evolve as a function of volume fraction.

All eight of the deterministic models (both KC and RGPZe) behave similarly and are clearly both significantly influenced by the variation of the mixed porosity, with the minimum in permeability controlled by and matching that of the mixed porosity. Moreover, the mixed porosity exerts a control on permeability leading to a variation in permeability of about two and a half orders of magnitude for 0 < χVf < 1.

However, the absolute value of permeability is also influenced by the type of model used to calculate the effective grain diameter. This is to be expected as the permeability calculated from both the KC and the RGPZe model is proportional to the square of the grain diameter (Eqs. 34) and (35). Consequently, higher permeabilities at any mixing ratio are provided by the arithmetic effective grain diameter model, moderate permeabilities for the geometric model and low permeabilities for the harmonic model. The cubic scaling model provides permeabilities similar to those provided by the harmonic model at χVf < 0.1 and similar to those provided by the geometric model for χVf > 0.6 following the behaviour observed in Fig. 9. The importance of knowing which of the effective grain size models to use is underlined by the recognition that there is a variation of up to about one order of magnitude according to which model is used.

6.4 Experimental Data

Figure 8 also includes one set of experimentally determined permeabilities, in order to provide some practical reference for all of these models. Qualitatively, the experimental data seem to be fitted well by a number of models. The RGPZe/geometric mean effective grain size model (solid green line in Fig. 8) seems to perform best, but other models, such as the Kozeny–Carman/Cubic combination, also appear to do well at some mixing fractions. We have calculated a statistic to judge which of the models best fits the experimental data. We use the sum of the squares of normalised differences between the experimental data and the theoretical prediction for each theoretical curve (SSND). The normalisation takes the form of dividing the difference between the predicted and experimental data by the experimental data and is described in more detail in Table 3. This normalisation is a modification of the common sum of squares statistic to ensure that all data are equally weighted even though the permeability varies by three orders of magnitude. In effect, the normalisation ensures that the effect of many good fits at small permeability values is not swamped by one bad value at high permeabilities.

Table 3 Example SSND calculations for fit testing for Mixtures 8 and 12

Table 3 contains an example of the SSND calculation for the fitting of models to the data for Mixtures 8 and 10 as shown in Figs. 6a and 8. The same data are shown in Fig. 10. Table 3 and Fig. 10 show that the RGPZe model with a geometric mean calculation for effective grain diameter provides the best overall fit to the experimental data (SSND = 0.864), with the KC model using cubic scaling being second best (SSND = 1.231). The worst fit was provided by the RGPZe with harmonic mean calculation of effective grain diameter (SSND = 62.01) and KC with harmonic mean used to calculate the effective grain diameter (SSND = 12.767) models. Similar SSND calculation tables for the other mixtures show similar results.

Fig. 10
figure 10

The square of the normalised difference \( {\text{SND}} = \left( {\left( {k_{\text{expt}} - k_{\text{model}} } \right)/k_{\text{expt}} } \right)^{2} \) as a function of χVf for each model. Low values indicate a better fit of each model to the experimental data for each mixture

When theoretical curves were compared with the other 11 binary sets measured in this work, it was found that the RGPZe model with a geometric mean calculation for effective grain diameter provided the best fit in all but three cases, and in all of these three cases, the RGPZe/geometric model came very close second. Those three exceptional cases were for the Mixtures 8 and 11, Mixtures 9 and 11 and Mixtures 9 and 10, the former of which the Kozeny–Carman/Harmonic model worked marginally better and the latter two of which the RGPZe/arithmetic model worked marginally better. In order to provide a more quantitative overall judgement on model performance, we have calculated the sums and arithmetic means of the SSNDs over all 12 mixtures for each model tested, which are ranked in Table 4, the best towards the top, having the smallest mean SSND. Although many of the models provide perfectly valid approximations to the data for some of the experimental measurements, the best models overall are, in order of decreasing efficacy, (1) RGPZe/geometric, (2) RGPZe/arithmetic, (3) KC/cubic and (4) KC/geometric.

Table 4 The arithmetic mean of SSND over all 12 mixtures studied in this work: the lower the value to better the fit

Figure 11 shows the measured permeabilities (the data given in Table 3) for all combinations of the bead and sand sizes shown in Table 1 as a function of volume fraction of fine beads or grains χVf. The uncertainty values on the experimental data vary with permeability and hence with flow rate as discussed in “Appendix”. The errors have been included in the figure by using the data given in “Appendix” to produce an overall percentage error as a function of permeability for 41 different permeabilities varying over five orders of magnitude. A simple linear curve can be fitted to these data with a coefficient of determination better than R2 = 0.98. The equation is given by

$$ e = 1.57 \times 10^{ - 4} k + 4.82, $$
(40)

where e is the total uncertainty (in  %) and k is the measured permeability in (D). This equation has been used to allocate a y-axis uncertainty to each point in Figs. 8 and 11. The x-axis uncertainty has been calculated to be constant at 1.024%, and this is applied directly to Figs. 8 and 11.

Fig. 11
figure 11figure 11

Measured permeability as a function of volume fraction of fine beads or grains from experimental measurements (symbols with calculated uncertainties). The solid lines represent the four different implementations of the RGPZe model, while the dashed lines represent the different implementations of KC model

Each part of Fig. 11 also includes theoretical curves for the model that provided the best fit, and in the case of the Ottawa sand and fine sand mixture, for the best fit theoretical curve and the RGPZe/geometric mean combination. In all cases, the measured permeabilities did not conform to the weighted interpolations between the end-member permeabilities

In summary, in almost all cases the permeability of the mixtures accords well with the theoretical curves that have been generated using the porosity equations developed in this work (Eq. (29)) together with the RGPZe equation, the geometric mean method for generating an effective grain diameter for the mixture and a simple mean value of end members for calculating the cementation exponent.

7 Conclusions

A study of binary mixing has been carried out involving the development of equations for mixed porosity and permeability from first principles. A total of 165 experimental determinations of porosity and permeability for a range of grain/bead diameters and volume fractions have been carried out, determining the accuracy of the measurements with precision.

The main conclusions are:

  1. 1.

    The porosity of a binary mixture can be characterised by two process: an interstitiation process where small grains fill interstices between the larger grains and a replacement process where large grains replace smaller grains and the porosity existing between them.

  2. 2.

    The equations describing these processes are independent of grain size, relative grain size, grain shape and the degree of packing. The effect of these four characteristics is exercised through their control on the end-member porosities, i.e. the porosities of the sample composed solely of the larger or smaller particles.

  3. 3.

    The porosity has a local minimum which is equal to the product of the porosities of each of the two size components, agreeing with the classical result.

  4. 4.

    The porosity minimum occurs at a grain fraction of small grains that depends only upon the end-member porosities.

  5. 5.

    Experimental measurements on 12 binary mixtures made from five sizes of glass beads and two grades of sand, each of which was tested in between 12 and 17 volume fractions of the fine particles, showed all measurements conformed to the theoretical equations.

  6. 6.

    Experimental measurements of permeability showed that conventional weighted mixing models cannot account for the permeability measurements of the binary mixtures.

  7. 7.

    Eight methods for calculating permeability were tested. All methods are dependent on grain size, porosity and the connectedness of the pore space.

  8. 8.

    The experimental measurements showed permeabilities that were clearly affected by the porosity minimum in the porosity mixing model.

  9. 9.

    The experimental measurements were described best by theoretical curves generated by using the exact form of the RGPZ equation, with the porosity from the porosity model developed in this work, together with a weighted geometrical mean grain diameter and using a simple arithmetic mean of end members to calculate the cementation exponent of the mixtures.

The model, as presented in this paper, assumes that grains touch at a single point to all other grains which surround it. There is no account taken of the compressibility, deformability, any degree of cementation or dissolution. Nor is there any account taken of the composition of the grains, clay swelling and wettability. Clearly, some of these secondary processes will reduce porosity, while others will increase it. One aspect of future work might be to modify the equations given above in order to take account of some of these processes.

It was the purpose of this paper to provide very-high-quality data that would either validate or refute the theoretical development that this paper describes. However, there are data within the scientific literature concerning geological and non-geological binary mixtures or near binary mixtures to which the equations developed in this work might be applied as future work.

This paper concerns itself entirely with theoretical developments and experimental validation. Other future work might include the implementation of 3D digital rock models (DRMs) to confirm the theoretical equations and check consistency with the experimental data. The advantage of such models is that they can more easily be modified to take account of secondary processes, say cementation, which might guide the modification of the basic equations to take account of each secondary process. Furthermore, the results of DRMs can be analysed to provide coordination numbers and other percolation information for the system.