Interactive comment on “ Percolation properties of 3-D multiscale pore networks : how connectivity controls soil filtration processes ”

The paper addresses the problem of straining efficiency in filtration processes. The model proposed exploits the multiscale nature of the pore size distribution, extending to the 3D case a method published by two of the Authors (Bird and Perrier 2009) for the 2D case. The Authors use (1) analytical renormalisation methods and (2) numerical methods, to study the CFS (critical filtration size) as a function of the structural and geometrical parameters of the network. The paper is relevant to Biogeosciences and the results are worth being published. However, the Authors should consider the following points: (1) In the Introduction (3001, line 18), the Authors state "we provide an improved algorithm for numerical simulation". The algorithm is never explained in the paper. On page


Introduction
Filtration of impure water by soils has been studied from many points of view and at many different scales.Filtration processes involve several biological, chemical, and physical mechanisms which are directly or indirectly linked to pore sizes and whose exhaustive study is far beyond the scope of the present paper.In the present paper we will not address all the components of filtration processes, but only the straining component."Straining is a matter of particle size alone, with the particle either larger or smaller than the pore through which it is attempting to pass" (AFS, 2009).Bradford et al. (2006Bradford et al. ( , 2009) ) report recent experimental evidence that indicates that straining may explain many of the reported limitations of former filtration theories that did not include the potential influence of physical factors such as soil pore size distribution and surface roughness and that the latter can play an important role in colloid deposition and microbe retention.
Our objective is to propose a new model of straining efficiency based on the knowledge of the multiscale nature of the pore size distribution observed in many soils, and on new developments on percolation properties in multiscale pore networks (Bird and Perrier, 2009).We will consider in this methodological paper any type of particles so-called "contaminant entities" that we would like to trap inside the soil (bacteria, viruses, fungi, any type of colloid, etc.).
A first approach to straining processes could be monoscale.Let us consider a simplified, water-filled porous medium, with pores of a given linear size r: entities of size bigger than r will not be allowed to enter the soil (surface filtration) whereas all entities of size smaller than r can potentially be dispersed by water inside the soil.If there exists a path for water to percolate through the network of pores of size r, and if we neglect clogging processes, the same paths will allow water flow and the convective transport of entities through the soil.Conversely if no path occurs, the porous medium will become not only a perfect filter for all entities (depth filtration) but also impermeable to water flow.
More realistic models of depth filtration (Price et al., 2009) account for the irregular 3-D geometry of the pore space, by introducing large pores and small throats.McGechan (2002) also introduces a broad range of pore sizes whose distribution is derived from hydrodynamical properties, in a mathematical model, and he concludes that large entities filtered by straining might be very rare.On the contrary, Bradford and Bettahar (2005) consider from experimental grounds "that straining may play a much more significant role than previously thought".
In real soils, visualisation techniques can give some information at observable scales.In their review, Keller and Auset (2007) report recent results obtained at the pore scale which show that some filtration processes are colloid-size dependent (the reduction in pathways for biocolloids as a function of their size lead to earlier breakthrough) and some other ones may be flow-velocity dependent (one observes that the largest entities are located in the central flow streamlines and the smallest at the interfaces) contrary to previous knowledge.Let us note that the velocity of fluid flow depends in turn on the pore size distribution and on the connectivity of the pore network.Keller and Auset (2007) conclude that the interception of biocolloids by the porous media are so far poorly understood at the microscopic scale despite the increasing power of observation techniques during the last decade.
In the present theoretical paper, we will neglect all dynamical properties of filtration processes to focus on the link between the size of filtered entities and the size of pores, using a simple percolation model in a new, multiscale approach.
It is well-known that the range of pore sizes occurring in most soils is very large, and fractal geometry has been widely used to quantify a pore size distribution by means of a power-law function with few parameters: e.g. the largest pore size and a fractal dimension in a mass fractal model (e.g.Rieu and Sposito, 1991;Perrier et al., 1996;Bird and Dexter, 1997).The range of living or inert "entities" which can be encountered in soil ecosystems is of the same order of magnitude for obvious geometrical reasons (Marilleau et al., 2008;Blanchart et al., 2009).By the way let us note that McGechan (2002) wrote: "The fractal approach may have the potential to indicate irregularities in pore shape which will restrict colloid movement, but such a procedure has not so far been attempted."Let us also note that the irregularity of the pore-solid boundary is well described by mass fractal models, where the actual "geometrical fractal object" is the fractal interface of dimension D, and the so-called "fractal pore size distribution" is a mere power-law distribution whose exponent depends on D (Perrier et al., 1999).
The knowledge of the pore size distribution alone (e.g. from 2-D image analysis) is not enough to predict soil hydraulic properties: water flow strongly depends on the connectivity of the 3-D network made of interconnected pores of different sizes (Perrier et al., 1995;Bird and Dexter, 1997).In the same way, the knowledge of soil depth filtration properties requires an estimation of available pathways inside soils (Price et al., 2009).
The issue can be considered as a double connectivity issue.Water will percolate through a soil sample if there exists a continuous path connecting one side to the opposite one.Entities will be filtered by the soil sample if there exists no continuous path of sufficiently large pores connecting their entrance location to the exit one.The critical size for filtration (CFS) will be defined as the size of the largest contaminant entities which will not be trapped in the soil because they will find a continuous path of pores of size CFS to convey them towards the exit.
One might address the issue in 2-D, as did Kaiser (1997) for directed percolation for clogging in a porous medium, but we have shown (Bird and Perrier, 2009) that realistic models of soil connectivity can be achieved only in 3-D.So we will work in 3-D, on the case study of a 3-D mass fractal.We will show here that previous results obtained in 2-D can be extended to the 3-D case using new renormalisation functions and we provide an improved algorithm for numerical simulation.
The first section will present theoretical arguments based on renormalisation functions to show qualitatively how filtration depends on the structural-geometrical parameters of the 3-D model.The second section will present the results from numerical experiments carried out on 3-D prefractal media to evaluate the CFS value as a function of structuralgeometrical parameters.The final section will discuss how far these first results open the way towards a new methodology to estimate soil filtration efficiency from the construction of soil structural models to be calibrated on available multiscale data.

Theory
Standard percolation theory in random monoscale 3-D networks has been established long ago and details can be found in textbooks (Stauffer and Aharony, 1992).
In a simple cubic site percolation lattice, sites are occupied with pores with a fixed probability p.The percolation threshold, which delimits connected from disconnected large networks, has been determined by numerical simulations up to the 7th decimal, (e.g., Tiggermann, 2001).
The existence of such a threshold can be proved using renormalisation theory (e.g., Lesne, 1998) and an easy estimation of its numerical value can be obtained as the fixed point of a renormalisation function.For a 3-D cubic site percolation network, the renormalisation function is (Turcotte, 1992): and the threshold estimate is given by the solution of f (p) = p, that is: Unfortunately, as mentioned in the introduction, the monoscale approach does not suit well the reality of soils.Bird and Perrier (2009) recently applied a renormalization approach to a multiscale fractal network, and their first results were analysed in the 2-D case.The generalization to the 3-D case is straightforward.As in the 2-D case, the new renormalization functions for the multiscale pore and solid networks become: where q is the scaling parameter of the mass fractal model (see next paragraph) and f (p) is given in the 3-D case by Eq. ( 2a).Repeated application of these functions through yields the change in probability of connection in each phase as more levels of structure, i, are added to the fractal (Fig. 1).
The classical random mass fractal model we adopt is well documented in the soil science literature (Bird and Dexter, 1997;Rappoldt and Crawford, 1999;Sukop et al., 2002).It is based on an n×n×n grid.We choose the simplest case, namely n = 2. Pores are assigned with probability q to each grid element.Those which remain unoccupied are then subdivided using a finer grid, and smaller pores assigned as before in a recursive way.The resulting structure can be characterized by its fractal dimension given by: and its porosity after i iterations is given by: Adding levels in a mass fractal model increases the porosity : if the number of levels i tends towards infinity, the porosity tends towards 100%.A realistic model involves an upper cutoff of scale associated to a maximum pore size r max as well as lower cutoff of scale r min , that is a finite level of iterations in a so-called prefractal model (Sukop et al., 2002).The probability to percolate for a given fractal model of parameter q is also an increasing function of the number i of iteration levels, but in a non-trivial and non-linear manner which has been depicted in previous studies (Bird and Perrier, 2009).
From a qualitative point of view, the 3-D model behaves in a quite similar way, but the connectivity of a real soil is better modelled and it has been shown that the 3-D approach is mandatory to account for plausible connectivity values for both the pore network enabling transport and the solid network enabling structure stability.
Figure 1 depicts graphically the iteration process given in Eq. ( 4) for both the solid and pore phases.Through this we may observe how the connectivity of the network evolves for different values of q (and hence D).The renormalisation functions f pore and f solid are represented respectively by blue and brown curves.For the pore phase we start with a probability of occupation p = q.We then pass vertically to the curve f pore (p), yielding p .We then cross to the diagonal line p = p, effectively swapping p for p.We may then repeat the above procedure, resulting in iteration of Eq. ( 4).Similarly for the solid phase we start with a probability of occupation p = 1 − q and use the curve f solid (p).Both iterative schemes exhibit fixed points where the function f intersects the diagonal p = p.For high values of q (Fig. 1a, q = 0.4) and low values of D, one can see that only 3 levels enable the pore network to percolate with a probability close to 1 (0.99998) whereas the probability to percolate for the solid phase drops towards zero (0.00096) with only 7 iterations.Such a structure clearly has no relevance in terms of modelling soil.For a critical value of q (q = 0.28011, see Fig. 1b), the diagonal becomes tangent to the f solid curve.This is called a tangent bifurcation and leads to a new fixed point to which iterates are attracted.The probability for the solid phase to be well connected is then high for an arbitrary number of fractal iterations.When q decreases from this first tangential value towards 0 (from Fig. 1c to f), the probability for the solid network to percolate always remains high and the soil structural model is very stable.
As regards "entities filtration" as well as "water infiltration" processes, we now focus on the connectivity of the 3-D pore network enabling water flow and particle transport.When q goes on decreasing (Fig. 1c to e), one passes through a critical transition when the diagonal becomes tangent to the f pore curve.Before reaching the tangent bifurcation (at q = 0.06917), one needs more and more levels for the pore network to percolate.
Once we have reached the tangent bifurcation, the probability of pore connection approaches a new fixed point which is too small to represent a viable model of soil, irrespective of the number of levels of pore structure included.Indeed, as in 2-D previous studies, we can achieve in this regime structures with porosities arbitrary close to 100% with arbitrary low connection probabilities.
More generally, the result of the theoretical analysis on the filtration issue is the following: Filtration is linked to percolation of the pore network.In the mass fractal model, the type of percolation behaviour of the pore network does not depend on the value of porosity but on the critical parameter q (associated with a given fractal dimension through Eq. (5a) www.biogeosciences.net/7/3177/2010/Biogeosciences, 7, 3177-3186, 2010  3) for the pore network (blue curve) and the solid network (brown curve).Figure 1a to 1f depict the results of the upscaling renormalization process for varying q, each step in the step functions plotted respectively in red or brown represents an iteration level in the fractal model using Eq. ( 4).
We will see in the following section how a critical CFS can be defined by the pore size associated to the minimum number of iterations enabling the probability to percolate to reach a value beyond p c = 0.31... (Eq.1), which is to the value enabling a random array of fractal cubes to percolate.
The present theory is mainly qualitative, because the wellknown discrepancy between the reference thresholds calculated by means of computer values and the values calculated from renormalisation (resp.Eqs. 1 and 2b) increases also with the number of embedded levels in the model.The next section will give the results of numerical experiments designed from the present theoretical analysis.

Percolation algorithm
A novel algorithm has been developed to better handle percolation in large 3-D fractal structures by accounting for the presence of large clusters of connected voxels.
A classical algorithm is based on the initialisation of a 3-D array of conducting or non conducting cells, and on the exhaustive search for connected conducting cells from one side of the array to the opposite side.Such an algorithm was used by Bird and Perrier (2009) for 2-D networks.Unfortunately this algorithm is severely limited by memory, as it is necessary to use large arrays to record the status of the network.
The new algorithm is more efficient in terms of memory allocation.It works in successive steps of an upscaling scheme based first on the detection of conducting clusters from a local search of neighbours around each conducting cell, then on the merging of the connected clusters.The whole array percolates when two cells from opposite sides belong to the same infinite cluster.This scheme enables a partition of the whole cubic array into n 3 subcubes.The upscaling cluster scheme can be applied first in any subcube, then the interior cells can be erased from the computer memory because only the boundary cells of each subcube are relevant when the whole array is reconstructed.One can work on the subcubes independently, enabling parallel computing.The last improvement of the algorithm then consists in generalizing the decomposition into subcubes, using a recursive, downscaling approach.This is well adapted to the fractal model because some subcubes can be at once identified as conducting clusters or non conducting clusters at each scale.Thus the new algorithm takes also advantage of the fractal hierarchical structure to avoid neighbours computer search and to spare some computer running time.
The results have been compared with the classical algorithm and the result of a simulation on any given array are identical (over 1000 tests).(By the way the estimation of the percolation threshold from a set of 1000 simulations on random cubic arrays is identical to the classical values found in textbooks).
In the future, because the hierarchical iterative space partition used by the novel algorithm enables the parallelisation of the code to be run on computer clusters, the new algorithm should become much superior as regards both memory handling and running time, and even in the random non fractal case.

Simulations
We carried out numerical simulations on prefractal cubic models of linear size L developed over 10 iterations, to represent a plausible range of scales as regards the modelling of many soil pore size distributions.Thus the largest pore size r max corresponds to i = 1 and the smallest r min to i = 10, which means r max = L n and r min = L n 10 .Taking L = 1 mm gives a range of pore sizes from millimeter to micron scale.
Table 1a gives the probability to percolate estimated from 100 runs of the percolation algorithm for n = 2 and values of q ranging from 0.2 to 0.01.For q = 0.2, one can check that the probability for the pore network to percolate reaches 1 from level i = 5, which matches very well the theoretical result given in Fig. 1c.For smaller values of q, one observes the expected discrepancy between numerical estimates given by the renormalisation approach and the direct computer simulations: The percolation regime occurring below the tangent bifurcation, (where even an infinite level of iterations will not allow percolation) appears only around q = 0.04... that is for smaller q values than the one predicted by the theory (around q c = 0.07).The results are nevertheless qualitatively similar.For small values of q (q = 0.03), the probability to percolate is very low (it equals only 0.06 at i = 10 iterations), whatever the total porosity value (the total porosity is here already 26%, Table 1b).One knows from previous theoretical arguments that further iterations would increase the porosity towards the 100% limit for infinite iterations whereas one does not expect percolation at all!As regards filtration, the principle is illustrated on Fig. 2. The visualization of a computer simulation is shown on an arbitrary statistical realization of a mass fractal model (n = 2, q = 0.1, thus D = 2.85 from Eq. 5a) developed over 10 iterations.The first 5 levels only are plotted in Fig. 2 for sake of visual resolution (Fig. 2a).From Table 1a, the probability to percolate reaches 1 at level i = 7, so we could predict that water will almost certainly go through this modelled soil as soon as i = 7. Percolation may happen at any level: Fig. 2b exhibits a continuous path of pores conducting water from one side to the other one at i = 5, where the probability of occurrence of such a path is 0.48 (Table 1a).Entities whose size is larger than the size of pores of level 3 have a probability to be trapped of 0.81, since the probability for the pore sub-network of level 3 to be connected equals 0.19 (Table 1a) and they are actually trapped on Fig. 2c.In a similar way, entities whose size is larger than the size of pores of level 4 have a probability to be trapped of 0.52 (1−0.48) and they are actually trapped on Fig. 2d.
The following question is: can we actually define a CFS value using such a statistical fractal model?The answer is the following: The CFS value is the size associated to the transition between a non percolating regime and a percolating regime, when the probability to percolate switches from a value below p c to a value above p c , where p c is given by Eq. ( 1).Why?A representative soil sample cannot be modelled by a simple (pre)fractal cube with a few pores of size r max .Thus we can design an improved model of a real soil by means of the juxtapostion of random fractal cubes, that is a random 3-dimensional array of fractal cubes.To actually build this array, the maximum pore size of each fractal cube has to be selected as the maximum pore size occurring in the real soil (by image analysis or from the air entry value in water retention data).The second parameter of the fractal model is its dimension D (calculated from the pore size distribution scaling).Each pre-fractal-sub-cube has a probability p f to percolate given in Table 1a.Then, from Eq. ( 1), the random array of fractal sub-cubes percolates when this probability p f is higher than p c .
A graphical solution is given by Fig. 3.For q = 0.05, we have 7 filtering (non percolating) levels, and 2 more levels which permit water flow through the pore network and transport of particles of size smaller than a CFS corresponding to level i = 8, (4 micron).For q = 0.2, we have two 1 filtering levels, and 8 more levels permitting water flow and transport of particles of size smaller than CFS corresponding to level i = 2 (0.25 mm) The higher q (and the smaller D for the same n = 2), the higher the connectivity, but the larger the CFS or the smaller filtering efficiency.
The model is discrete.The CFS values are found among the sizes associated with each iteration level i.They can be visually estimated from Fig. 3 or numerically from two successive gray cells in any row of Table 1.Conversely, the actual critical q c values where the transition between the nonpercolating regimes occurs are estimated by a simple linear interpolation in the q-interval defined by two successive yellow cells in any column of Table 1 or by inverse simulations by means of numerical trials and errors on q values.The resulting q c values are given in Table 2, as well as associated D c values from Eq. (5a).
The associated changes in CFS as a function of D are plotted in Fig. 4, which shows the good qualitative agreement between renormalization theory and numerical experiments.1a simulated results) for varying values of parameter q in the fractal model (and associated values of fractal dimension D and porosity calculated in Table 1b).The straight line p f = p c (Eq. 1) delimits the transition between the percolating and non percolating regimes for a random array of fractal cubes.
Table 2. Critical q c estimated by a linear interpolation in the [q1, q2] interval defined by two successive yellow cells in each column of Table 1 and probabilities pf q 1 and pf q 2 .D c is calculated from qc using Eq.(5a).

Discussion and conclusion
The results of the present study are both theoretical and methodological.
From the theoretical point of view, we have shown that the multiscale percolation approach developed by Bird and Perrier (2009) extends to the 3-D case in a straightforward way, and that a 3-D mass fractal model can allow for both percolation of the pore network to enable fluid flow and percolation of the solid phase to ensure stability.
From the methodological point of view, we have shown that we can define a critical filtration size CFS from two types of information: information on the pore size distribution and information on the connectivity of the pore network.Information on the pore size distribution is easily available from image analysis (2-D representative soil sections are enough) see Bartoli et al. (1999) for how to measure the pore distribution and also Perrier et al. (2006) for how to compute the mass fractal dimension for an 2-D image of the pore-solid interface.Information on the connectivity of the pore network is more difficult to obtain in a direct way.3-D image analysis (Gryse de et al., 2005) offers some potential but the present state of available tools and their resolution excludes access to the whole range of pore scales.One could also try to work, in an indirect way, on a pore network calibrated by inverse modelling of hydrological properties depending on soil connectivity (e.g.Johnson et al., 2003;Price et al., 2009) but again neglecting the whole range of scales involved in water flow and particle transport.The present approach is based on an alternative approach, that is the construction of a very simple 3-D model developed with a multiscale approach.The methodology is based here on a widely used model of soil geometry, the mass fractal model, and it can be summarized by the four following steps: 1. Using image analysis to determine the pore size distribution and r max .3. Building a tridimensional multiscale pre-fractal model of dimension D to account for pore network connectivity at various scale.

Estimating CFS.
As regards potential use on experimental data, let us conclude here with general considerations: Bradford et al. (2003) mentioned that "Few studies have examined the influence of soil pore size distribution characteristics on colloid straining. . ." and that "they were unaware of any models that explicitly account for straining".They first estimated the percentage of pore space smaller than a critical straining pore size from capillary pressure curves which depend indirectly on connectivity (Bradford and Bettahar, 2005).We propose here a new methodology which explicitly accounts for connectivity by means of a geometrical model of a multiscale pore network, and which could be easily tested using complementary data using different sizes of contaminant entities in soils with a broad pore size distribution.The distribution of the solid particles sizes has also been found to have a high impact on straining (Xu and Saiers, 2009;Diaz et al., 2010).Diaz et al. (2010) showed that, according to experimental observations, "straining was shown to contribute highly to bacterial retention in all the soils that they tested, in particular in the soils with a broader grain size distribution and more irregular shape".Future modelling work will address the scaling relationships between pore and solid size distributions.The popular mass fractal model is but one example of a multiscale porous structure.Others exist, including the PSF model (Perrier et al., 1999;Bird et al., 2000) which accounts for broad ranges in size of both soil pores and soil particles, as well as for the irregularity of the solid-pore interface.Again models such as these can be substituted into the above methodology to provide alternative approaches to estimation of straining efficiency.We consider this work opens the way towards a new methodology to estimate straining efficiency in soils using complex networks calibrated on available multiscale data.

Fig. 1 .
Figure 1 2 3 (a) The pores are coloured white and the solids brown.Water coloured blue was injected from the hidden vertical side, and connected paths occurred when some cells are coloured blue on the opposite front side.(b)Visualisation of the cluster of connected blue "water-filled" interior cells of the 5th size class which defines the CFS of this example.(c)Visualisation of the cluster of green interior cells of the 3rd size class, it is not connected to the exit and contaminant entities with a larger size would be trapped.(d)Visualisation of the cluster of connected green interior cells of the 4th size class, which again would retain trapped contaminant entities.

Fig. 2 .
Fig. 2. A 3-D illustration for a random realisation of a fractal cube defined by n = 2, q = 0.1 and i = 10 iterations.Only the first 5 levels and associated 5 pore sizes are drawn.No pores of the 1st level occurred due to low q value.

Fig. 3 .
Fig. 3. Numerical experiments: Probability pf for the pore network in the fractal cube to percolate as a function of iteration level i (plotting Table1asimulated results) for varying values of parameter q in the fractal model (and associated values of fractal dimension D and porosity calculated in Table1b).The straight line p f = p c (Eq. 1) delimits the transition between the percolating and non percolating regimes for a random array of fractal cubes.