Flux-conserving directed percolation

We discuss a model for directed percolation in which the flux of material along each bond is a dynamical variable. The model includes a physically significant limiting case where the total flux of material is conserved. We show that the distribution of fluxes is asymptotic to a power law at small fluxes. We give an implicit equation for the exponent, in terms of probabilities characterising site occupations. In one dimension the site occupations are exactly independent, and the model is exactly solvable. In two dimensions, the independent-occupation assumption gives a good approximation. We explore the relationship between this model and traditional models for directed percolation.


Introduction
Percolation problems were introduced by Broadbent and Hammersley in 1957 [1].Their paper motivated the study of percolation [2] by a discussion of fluids passing through a disordered medium, such as water penetrating through limestone.Forced flow of a liquid through a porous medium is central to many interesting and technologically important processes involving elution from, or absorption by, random media, such as leaching of salts from soil, extraction of oil or gas from reservoirs, brewing coffee, or the operation of chromatography columns, which stimulated the study of directed percolation.There is a vast literature treating the standard models of directed percolation, reviewed in [3], and [4,5,6,7,8,10,9,11,12,13,14,15,16] are indicative of the breadth of different approaches to directed percolation, and of its wide range of applications.The standard percolation models (such as bond percolation on a lattice) do not take into account the mass conservation of a flowing liquid.In this paper we consider a generalisation of directed percolation, which includes the flux in a bond as a dynamical variable.If the fluxes are ignored, and only the occupancy of bonds is considered, then the standard directed percolation model occurs as a particular special case.Our generalised model includes a flux-conserving case which can describe the forced flow of a liquid through a random medium, and we shall consider this in some detail.We demonstrate that aspects of the subset of models which represents elution by a flux-conserving fluid are exactly solvable in one dimension, extending some results obtained (in another context) in [17].Our generalised model has some similarity to the Scheidegger model for the distribution of river catchment basins [18,19], which will be discussed in the conclusions.Another, more distantly related, class of models which quantify directed transport in random media is described in [20,21].
The percolation model introduced by Broadbent and Hammersley uses the idea of 'wetted' bonds.Here we extend this binary notion of wetting by modelling the flux ϕ of liquid through each bond and the probability distribution of such fluxes.We argue that, for a quite general class of flux-conserving process which involve a forced flow through a disordered medium, the distribution of the fluxes has some universal characteristics.After penetrating a sufficient distance into the network, the probability density P (ϕ) for the flux ϕ in a channel may approach a stationary distribution, with a power-law form at small values of ϕ: (1) The existence of power laws is usually associated with critical phenomena, however the power-law distribution described by equation ( 1) is a robust feature, which does not depend upon tuning the model to criticality.We argue that the mechanism determining the power law is very general, and that power-law distributions of small fluxes are a robust feature of the model.Section 2 defines our model, in both one and two dimensions.Our model includes both a 'skeleton' of wetted bonds, and the flux of material, ϕ, carried in each wetted bond.The skeleton is defined by three probabilities: a wetted channel continues with probability p 1 , splits into K channels with probability p 2 , and terminates with probability p 0 .The flux carried by a bond which splits is distributed so that a fraction r k flows into each branch, with r 1 + . . .+ r K = 1, where r k ∈ [0, 1].We emphasise cases where there are only K = 2 branches, and where the two branches carry fixed fractions of the total flux, denoted by r and 1−r.Accounting for the relation p 0 +p 1 +p 2 = 1, the model then has three independent parameters, p 0 , p 2 and r.The flux-conserving case is the p 0 = 0 subspace.The model which was considered in [17] is the one-dimensional, flux-conserving case.
Section 3 considers the distribution of fluxes in the mass-conserving case in some detail.It is shown that, for small values of the flux, the probability distribution function (PDF) of the flux is asymptotic to a power law described by equation (1).We obtain an exact equation for the exponent α of this power law, in terms of some probabilities, P j which characterise the occupation of lattice sites.
The bond skeleton is characterised by the probability f that a given bond is occupied.In section 4 this is calculated under the assumption that the occupation of sites is statistically independent, and we find very good agreement between theory and numerical experiment.Section 5 estimates the probabilities P j which define the equation for α using the same approach, and we compare empirically determined values for the exponent α with those obtained from both numerical and theoretical estimates of the P j .
In section 6 we consider the extent to which the independent-occupancy approximation is exact.We demonstrate two results which indicate that it is exact in the one-dimensional version of the model.We find that, in the two-dimensional case, this is a very good approximation, but not exact.
Section 7 discusses the relationship between our system and the standard model for directed percolation.By varying the parameter p 0 we can induce a percolation transition in our model, which can be regarded as a consequence of large voids appearing between the active bonds.In sub-section 7.1 we investigate where the transition lies in the parameter space of our model.Our model can be thought of in terms of a combination of 'bond' and 'site' deletion processes, and section 7.1 also discuss the relationship between our system and a model for mixed site and bond percolation, discussed in [16].Sub-section 7.2 presents some numerical evidence that the the critical exponents describing the structure of the skeleton are the same as for the usual directed percolation model.In section 7.3 we present some results on the distribution of void sizes.
Finally section 8 summarises and discusses the implications of our results, including a model for the slow elution of material by percolation.

Definition of the model
We define discrete dynamical processes, with an iteration number, j.Incrementing j can be thought of as advancing time, but in the physical contexts described in the Introduction, increasing j represents moving downstream in the forced flow.Because The model is defined on a bi-partite lattice.Only even sites (red) are occupied at even-numbered iterations, only odd sites (blue) during odd iterations.At each iteration every occupied site is either removed (probability p 0 ), moves to one of its nearest neighbours (probability p 1 ), or else splits into K = 2 daughters (probability p 2 ).The flux ϕ carried by an occupied site moves with it, unless the site is removed, in which case the flux disappears.When two fluxes are moved to the same site, they are added.
the iteration index can be interpreted as a discrete time, our one or two dimensional systems are comparable to standard models for directed percolation in 1 + 1 dimensions or 2 + 1 dimensions, respectively.

One-dimensional model
We consider a bi-partite lattice.At even iteration number j, only even sites are occupied.The dynamics always moves an occupied site by one unit left or right, so that for odd iteration number, only the odd sites may be occupied.Every occupied site, index i, is associated with a flux, ϕ i (j), at iteration j.
First consider the dynamics of the site occupations.At each iteration, an occupied site is annihilated with probability p 0 , or moves either right or left with equal probabilities p 1 /2, or else it splits into K = 2 branches with probability p 2 .
The corresponding dynamics of the fluxes is as follows.If a site is annihilated its flux disappears.When transitions bring two occupied sites to the same position, their fluxes combine, and if the site branches, then its flux is divided.So values of the flux are changed by two possible processes: coalescence : splitting : where either (r + , r − ) = (r, 1 − r) or (r + , r − ) = (1 − r, r), both cases with probability equal to one-half.The model is illustrated schematically in figure 1.More generally, we can make r ∈ [0, 1] a random variable, chosen independently for every bifurcation, with a PDF which is symmetric about 1 2 .The model only describes a volume-preserving flow when the stopping probability p 0 is equal to zero.This parameter is included for two reasons.The parameter p 0 is included primarily so that our system encompasses the standard model of directed percolation.In addition, the case where p 0 > 0 can model situations where the solvent disappears, for example by evaporation.
The fluxes at each wetted site can increase (due to coalescence) or decrease (due to splitting), and some sites become unoccupied.The long-time behaviour of the model is characterised by the probability, f , that sites are occupied and the distribution of the non-zero values of ϕ.A nonzero value of p 0 implies a loss of conservation, which significantly changes the nature of the problem.When p 0 is large enough, the probability that a site is occupied at long times becomes zero, and as p 0 decreases, one observes a transition when the probability of occupation becomes strictly positive.

Two-dimensional model
Consider a model defined on an 2N × 2N square lattice, with sides identified to make a toroidal topology.Every site (i, j) carries a weight ϕ ij .Initially only sites with i + j even are occupied, with unit weight, so that our model is defined on a bi-partite lattice.The value of N is assumed to be large.
The lattice configuration is then evolved in discrete timesteps.At each site, we choose (with equal probabilities) a move to one of the four nearest neighbour sites.With probability p 1 , all of the flux on (i, j) is moved to the selected neighbour.Alternatively, with probability p 2 , branching onto K nearest neighbours occurs.Note that these moves preserve the bi-partite property, so that at iteration k, only sites with i + j having the same parity as k may be occupied.
We investigated two different versions of this model, which we term the two-branch and four-branch models.In the two-branch model, we have K = 2, and a fraction 1 − r of the weight at (i, j) is moved to the selected neighbour, and a fraction r is moved in the reciprocal direction.All of the random choices are independent.The model has three parameters of interest, p 2 , p 0 and r.
In the four-branch model, there is branching to all of the nearest neighbour sites, so that K = 4.The fractional weights for each branch are defined by four numbers, {r 1 , r 2 , r 3 , r 4 }, with r 1 + r 2 + r 3 + r 4 = 1.

Power-law distribution of fluxes in the mass-conserving case
In the following, Sections 3, 4, 5 and 6, we consider exclusively the mass conserving case, with p 0 = 0. To simplify the notation, we will set, in these sections, p = p 2 , so that p 1 = (1 − p).
We argue here that the distribution of fluxes P (ϕ) has a power-law behaviour in the limit as ϕ → 0 (with the exponent in (1) satisfying α < 1, so that the distribution is normalisable).Branching of a channel reduces the flux due to multiplying by a random factor r k < 1, which we assume to have a known PDF.Because fluxes are added when coalescence of channels occurs, this process increases the flux.We are interested in the distribution of very small values of the flux ϕ.In this case splitting and coalescence have very different effects.In the case where a channel carries a very small flux, coalescence with another channel will produce a much larger flux (with a value which is typically Figure 2: In two dimensions, when the flux at an occupied site splits there may be up to four branches.We consider two cases in some detail.In the two-branch version of the model there are two branches, which go in opposite directions, either horizontal or vertical, with equal probability.In the four-branch model, where branching events reach all nearest neighbours, with the locations of the weights r 1 , r 2 , r 3 , r 4 randomly assigned. comparable to the mean flux).Almost all coalescence events will, therefore, remove a very small value of ϕ, whereas splitting events just reduce its value.The distribution of very small values of ϕ is, therefore, the result of a competition between two process: the small values of ϕ continue to decrease due to splitting, but they are annihilated by coalescences.
In order to explain why a power-law distribution of the flux is expected, we start by making a change of variables.Instead of considering ϕ, we consider the probability density function (PDF) of a logarithmic variable, ψ ≡ ln ϕ.Consider the dynamics of the variable ψ (regarding increasing ψ as a displacement to the right).With every bifurcation of a channel, the points representing the values of ψ are split into K new points, and each one is displaced by ln r k .When two channels coalesce, the two values of ψ are replaced by ψ = ln [exp(ψ 1 ) + exp(ψ 2 )].In the following we shall assume that the PDF P (ϕ) is bounded so that the probability of ϕ being less than ϕ 0 approaches zero as ϕ 0 → 0. This is consistent with the distribution (1) provided α < 1.Under this assumption, in the limit as ϕ → 0, most coalescences occur with channels carrying a much larger flux.As a consequence, coalescence of a channel with a small flux, ψ, is replaced by a value close to that which characterises a typical channel.This picture implies that the variable ψ drifts to the left with each bifurcation, but, in the case of small fluxes, coalescence almost inevitably causes a jump back to a position close to the origin.
Because the equations defining the dynamics of ψ become independent of the value of ψ in the limit as ψ → −∞, the PDF of ψ should reflect this translational symmetry.In the limit as ψ → −∞, the PDF of ψ should be asymptotic to an eigenfunction of the translation operator.Because the exponential function is an eigenfunction of a translation operator, we expect that the PDF of ψ has the form P ψ (ψ) ∼ exp(λψ) . (3) Figure 3: Results of numerical simulations demonstrating that the distribution of flux, P (ϕ), is asymptotic to a power law at small ϕ.We show results for the two-dimensional model.Upper row: two-branch system: each panel illustrates different values of p, with r = 0.01 (left panel) and r = 0.5 (right panel).Lower row: four-branch system, with r = 0.025 (left panel) and r = 0.25 (right panel).
Note that we must have λ > 0 to have a normalisable distribution if this law holds as ψ → −∞.The corresponding distribution of ϕ is then a power law of the form This is a very general argument indicating that the steady-state distribution of fluxes approaches a power law as we go deeper into the percolation medium, but it does not yield a prediction of the exponent α. Figure 3 illustrates numerical simulations of the distribution P (ϕ) for the two-dimensional model, demonstrating that it is indeed asymptotic to a power law at small values of ϕ.
To calculate α we shall determine a master equation for the PDF of the variable ψ, valid in the limit as ψ → −∞.Because coalescence almost inevitably results in the value of ψ making a large jump to the right, we must consider the fate of sites which are occupied and which have not experienced a coalescence for a large number of iterations.The latter requirement is imposed because it is only those sites which have very small values of ϕ.For this subset of sites we introduce the following probabilities: P 1 ≡ probability of moving without branching or coalescence P 2 ≡ probability of moving without branching and undergoing coalescence P 3+k ≡ probability of branching K ways, with k branches undergoing coalescence (5) Note that P 1 + P 2 = 1 − p and P 3 + . . .+ P 3+K = p.In practice, when we estimate the P k from numerical simulations, we only accumulate statistics for those sites which have not undergone coalescence in the preceding N thr iterations.We take a sufficiently large value for this threshold such that the P k are insensitive to the value N thr .
To describe the asymptotic form of the probability of very small fluxes, we write down an equation for the PDF of ψ at iteration j + 1, which is valid in the limit as ψ → −∞.Note that, because events involving coalescence almost always induce a large increase of the flux, they do not contribute to this balance equation for very small values of ϕ.It follows that it is only events which do not involve coalescence at iteration j which contribute to P ψ (ψ, j + 1), when ψ → −∞: these events are moving without coalescence, leaving ψ unchanged (probability P 1 ), or splitting events where some of the daughters escape coalescence (probabilities P 3 , . . ., P 2+K ).Taking this remark into account, if P ψ (ψ, j) is the PDF of ψ at iteration j, then with Seeking a solution of the form (3) which is independent of j gives an exact equation for the exponent α: In the cases where K = 2 and where the splitting ratios are (r, 1 − r)), (this includes the one-dimensional model and the two-branch model in two dimensions) equation ( 8) simplifies to with We also simulated a model with K = 4: when there is a four-way split, we set a weight factor of r for two of the sites (chosen at random), and a factor 1 − r for the other two sites.In this case, α satisfies equation (9) with (four − way splitting) .(11)

Occupation probabilities in the mass-conserving case
Here we present calculations for the occupation probability f in the mass-conserving case, p 0 = 0, assuming that the occupation of sites is statistically independent of their neighbours.For the sake of simplicity, we set p = p 2 , so p 1 = 1 − p.We limit the discussion to the two-dimensional models, because the one-dimensional case was treated in [17], where it was shown that the occupation probability is where the sub-index 1 indicates the one-dimensional model.

Occupation probability: two-dimensional case with double branching
We determine here the probability of occupation in the two-dimensional problem, first with K = 2 (two-branch model), f 2,2 (p), where the sub-indices indicate the twodimensional model with two branches and we recall that p = p 2 .Assuming that sites are randomly occupied with probability f 2,2 , we estimate the probability P empty that a site will be empty at the next iteration.Note that the probability that one of the four nearby sites makes a transition to reach this site is The probability of the site remaining empty is This leads to a cubic equation for f 2,2 : for the two-branch model with prediction from independent-occupancy approximation, equations ( 17) and ( 18).(b) Shows the fractional error, plotted against ln p. (We defined the fractional error by ∆f where f th and f num are, respectively, the theoretical and numerically determined values of f .) The cubic equation arising from (15) can, in fact, be solved by the method of Cardano.Within the interval p ∈ (0, 1), we have only one real root.By this method, the dependence of f 2,2 on p is given by: where The maximum of f 2,2 (p) is f 2,2 ≈ 0.9126 at p = 1. Figure 4(a) compares this prediction of the filling probability with the result of numerical simulation.The agreement is very good, but not perfect.In particular, we find that the fractional error is quite large at small values of p, as illustrated in panel 4(b).We were not able to determine whether the fractional error eventually approaches zero as p → 0.

Occupation probability: two-dimensional case with fourfold branching
The occupation probability for the four-branch model in two dimensions will be denoted f 2,4 .For this model the transition probability is Then, the cubic equation for f 2,4 is: Shows the fractional error (defined in the same way as for figure 4), plotted against ln p.
Making a comparison with equation ( 15), we see that f 2,4 (p) = f 2,2 (3p), so that the occupation probability in this case is and, when p ≪ 1, we have f 2,4 ∼ 8p.In this case, the maximum value is f 2,4 = 1 at p = 1. Figure 5 compares this prediction of the filling probability with the result of numerical simulation.Again, while the agreement is very good throughout most of the range of p, there is a substantial fraction error for small values of p.

Estimates of transition probabilities in the mass-conserving case
In section 3 we presented an exact equation, (8), determining the exponent α in terms of a set of probabilities P k , defined by equation (5).In this section we consider these probabilities, comparing numerical simulations with theoretical estimates, where these are available.We consider different models in turn (including the one-dimensional model, because these probabilities were not given in [17]).In all cases, the theory uses the assumption that the sites are independently occupied with probability f , as estimated in section 4. In section 6 we shall argue that this independent-occupation assumption is exact in the one-dimensional case.Accordingly we propose that the formulae for the P k are exact in one dimension.
In the case of the two-dimensional model with two-way splitting, we are able to estimate the P k analytically using the independent-occupation model.We find close agreement with values derived from numerical simulations.In the two-dimensional model with four-way splitting, we are limited to giving values of P k derived from simulations.

One-dimensional model
A given trail can evolve in several ways at each iteration, with probabilities defined by equation (5).If sites are occupied with probability f 1 (p), we find, using Eq. ( 12), that there is a transition probability for a trail coalescing with one or other of its two neighbouring sites, given by (the term (1 − p)/2 comes from the case where the neighbouring site does not divide and moves in the direction that creates a collision, and p comes from the case where the neighbouring trail divides).
The probability P 1 arises from the case where a trail does not divide, and does not collide: Similarly In the case where the trail divides, there are two independent chances for the trail to be annihilated, so the probability for both daughter trails to end is Similarly the probability for both daughter trails to survive is And because there are two ways in which one daughter trail can continue Using equations ( 23) to (27) in equation ( 9), we find that for the one-dimensional fluxconserving model, α is a solution of in accord with Equation (2) of [17].

Two-dimensional model with two branches
Next consider estimates for the transition probabilities for the two-dimensional model with branching into two opposite directions.Again, we use the assumption that the sites are independently occupied, with probability f 2,2 , as approximated by equations ( 17), (18).First note that a site moves to one of its nearest neighbours with a transition probability, given by equation ( 13), namely A site remains un-branched with probability (1 − p), and un-combined if there is no other transition into the final site from any of its three other neighbours.Hence The value of P 2 is then determined by noting that 1 − p = P 1 + P 2 : If the trajectory branches (with probability p), it avoids collision if none of the three nearest neighbours of each of the two new sites make a transition which lands there.Hence Similarly, P 5 is the probability of an event where neither of the two branches avoids coalescence with at least on one of its three nearest neighbours: To determine P 4 we can use P 3 + P 4 + P 5 = p to obtain Figure 6 compares these theoretical estimates for the P k with numerical simulations.The numerical simulations include only sites which had not experienced coalescence in the preceding N thr = 8 iterations (the parameter N thr was introduced in the paragraph below equation ( 5)).These simulations show that the true P k values differ slightly from the theoretical expressions, equations ( 29)-(33), which are plotted in figure 6.The small difference between theory and simulation becomes negligible as p → 0. We used the theoretical expression for f (p), equations ( 17), (18), when evaluating equations (29)-(33) for figure 6. Plotting the P k for simulations with N thr = 0 yields curves which are barely distinguishable from equations (29)-(33), indicating that the small discrepancy is due to the theory neglecting the requirement to exclude sites which have undergone recent collisions, rather than the error in equations ( 17), (18).
It is impractical to impose very large values of N thr because, as p → 1, very few sites satisfy the requirement to have undergone no coalescences in the last N thr iterations: for N thr = 8, we found that the probability of an occupied site satisfying this criterion falls from 0.64 at p = 0.05 to 0.0013 at p = 0.75.Simulations with N thr = 4 gave points which appear coincident with those for N thr = 8 when included in figure 6  Figure 7 provides direct tests of the theoretical prediction for the exponent α.We estimated α from simulations with a wide range of values of r and p.In the left panel of Fig. 7, we compare the theoretical values of α obtained from equation ( 9) against numerically estimates, obtained by directly simulating the model: here we used the values of P k obtained from simulations with N thr = 8.In the right panel we collapse all of the data points onto two plots of F (p) ≡ (1 − P 1 )/Q, one obtained using equations (29) to (33), the other using values of P k obtained from simulations with N thr = 8.The former shows small but significant deviations at larger values of p.The small dispersion between the symbols gives an indication of the accuracy of our determination of the exponents α.
We remark that, if the errors in the approximations underlying equations ( 29)-( 33) and ( 17), (18) are negligible as p → 0, we can determine the limiting value of F (p) as p → 0. Using these expressions in equation (10)  (34)

Two-dimensional model with four branches
In the case of the four-branch model, a theoretical calculation of the probabilities P k is considerably more difficult, although we can obtain formulae for P 1 and P 2 which are analogous to those obtained in the two-branch case, and find  The calculation of the other P k is complicated, because once a site has branched from one site to its four neighbours, one has to consider transitions from the eight sites adjacent to the four newly occupied positions.Four of them can only reach one of the four new positions, but four of them could possibly reach two positions.We can, however, assert that P 3 = p + O(p 2 ) and conclude that lim p→0 Numerical investigation of the probabilities P k for the four-way splitting model is difficult because, except when p is small, the proportion of sites which do not undergo coalescence events is very small.Accordingly, we confined our numerical investigations to cases where p < 0.5.Figures 8 and 9, illustrating investigations of the four-branch model, are similar to figures 6 and 7, but do not include theoretical predictions of F (p) ≡ (1−P 1 )/Q.

Tests of exactness in the mass-conserving case
In the Introduction, we mentioned that, in one dimension, the case where p 0 = 0 appears to be exactly solvable.Here, we argue that the steady-state probability for occupying N consecutive sites at step j can be written as a product of independent probabilities at different sites.That is, if s i is the occupation of the i th site, then we postulate the joint probability to be the product:   Figure 9: Testing the determination of α for the two-dimensional four-branch model, using (8), and the P k derived from simulations with N thr = 6.Left panel: there is satisfactory agreement between the empirical values of α and the values obtained from (8).Right panel: the data points onto a plot of F (p) ≡ (1 − P 1 )/Q, using values of P k derived from simulations (with N thr = 6), compared with the values of F (p) obtained from simulations which included all sites (i.e.setting N thr = 0): there is a significant discrepancy as p → 1.

A X B
Figure 10: At iteration j + 1, sites A an B are influenced by their nearest neighbours (illustrated here for the two-dimensional model) at iteration j.Correlations may result from the fact that both A and B are influenced by occupation of the 'key site', X.
where f 1 is the probability of occupation of a single site, given by ( 12).This result can be demonstrated by assuming that (37) holds at iteration j, and testing whether the joint probability given by Eq. 37 remains unchanged by the dynamics.Because sites which are not adjacent to each other are not influenced by common sites at the previous iteration, non-adjacent pairs are obviously independent.In subsection 6.1 we investigate the joint probability P 2 (a, b) for two adjacent sites, and show that this factorises if we make an assumption about the relationship between the occupation probability, f , and the splitting probability, p.This relation need not necessarily be the same, as the function f 1 (p) given by ( 12), and we shall distinguish it by denoting this function by f1 (p).The same approach is also used to determine functions f2,2 (p) and f2,4 (p) which would ensure that P 2 (a, b) = P 1 (a)P 1 (b) for the two-dimensional models.We show that f1 (p) coincides with f 1 (p), implying that the one-dimensional model is exactly solvable, but that this does not hold for the two-dimensional cases.
We are also able to give an inductive demonstration of a more general result concerning the dynamics of the one-dimensional model: in section 6.2 we show that the boundary between an occupied and an unoccupied region fluctuates diffusively, while the occupation probabilities within the occupied region remain statistically independent.

Condition for factorisation
Consider two adjacent sites at step j + 1.These were influenced by the configuration at step j through their nearest neighbours.There is only one site, which will be referred to as the 'key site', which can influence both of the sites at step j + 1 (see figure 10).So if we are seeking to establish whether sites are independent, this one site should receive special attention.This observation is also true in higher dimensions.
Consider the influence of the key site, X, upon its two nearest neighbours, A and B. A 'wetted bond' connection may (probability P1 ), or may not, (probability P0 ), be made from the key site X to site A. Let a = 0 or a = 1 indicate whether site A is (respectively) empty or occupied.Also let P (a|0) be the probability that A becomes occupied at step j+1, given that no connection has been made to A from X, and P (a|1) be the probability that A is occupied if a connection is made from site X.Clearly, P (a|1) = a, because A is definitely occupied if the connection is made, so that P (1|1) = 1 and P (0|1) = 0.With these definitions, we can write an expression for the probability of A being occupied at step j + 1.The probability P 1 (a) of A being occupied has a contribution from a term where no connection is made from X to A (with probability P0 ), multiplied by the probability for the independent event in which site A achieves occupancy a by connections from its other neighbouring sites.Adding another term, representing events where X does connect to A, we have We shall need expressions for P0 and P1 .These depend upon which version of the model we consider.
One-dimensional model In the one-dimensional case, we assume that sites are occupied with probability f at step j.So the probability of X being occupied and splitting is f p, and of X occupied and moving to the left without splitting is f Now consider the joint occupation probability of sites A and B. Going from iteration j to j +1, there is a probability P00 that there is no transfer of occupation from X to either A or B. The probability that X transfers to B and not A is P01 , and the probability to transfer to both sites is P11 .In the one-dimensional case, these probabilities are The joint probability at step j + 1 for sites A, B is From (42), the occupations remain independent at step j + 1 if all three of the following conditions are satisfied: In the one-dimensional case, these conditions are These three equations all imply the same relationship between f and p: which is the same relationship between f and p as arises independently from the calculation of the filling fraction in the independent-site approximation, equation (12).Note that this calculation didn't require P (a|0) or P (b|1), only the probabilities P0 , P00 , P01 .It is, therefore, easily extended to the two-dimensional models that we considered.

Two-dimensional model
For the two-branch model In this case, equations (44) are satisfied by f2,2 (p) = 8p (1 + p) 2 . (47) Similarly, for the four-branch model In this case, equations (44) are satisfied by f2,4 (p) = 16p (3p + 1) 2 . (49) In contrast with the one-dimensional model, in the two-dimensional cases the functions f (p), which satisfy the factorisation condition, do not agree with the functions f (p) which describe the occupation probability of the sites discussed in section 4. We conclude that the site occupations are not independent in the two-dimensional models.

One-dimensional case: finite intervals
We have shown that in the one-dimensional case, the independent-occupancy distribution is self-reproducing under iteration of the model in the mass-conserving case, p 0 = 0. We can also establish a stronger result on exactness of solutions of the one-dimensional case.Suppose that we know that at iteration j, the occupied region is bounded: it is known that sites n − and n + are occupied, and that everything to the left of n − and to the right of n + is empty.All the sites with a parity different from j are empty.Let us assume that all of the sites n with the same parity as j satisfying n − < n < n + are occupied with independently with probability f 1 , given by (12).What can be said about the joint distribution of occupancy at the next (j + 1) iteration?
The end values n − and n + both change by ±1, independently.At each end the boundary expands with probability (1 + p)/2 or contracts with probability (1 − p)/2.Let us assume that the shifts of the boundary have been determined, and consider the joint distribution of occupation of the other sites, conditional upon the new positions of the ends of the occupied interval.
We shall argue that if, at iteration j, the interior sites are independently occupied with probability f 1 , then, they remain independently occupied with f 1 at iteration j + 1.Because of the short range of influence of occupation probabilities at the next iteration, we can treat the two ends independently, and consider only a short interval in the vicinity of one end.
We consider a sequence of four sites at iteration j, having the same parity as j.The first two of these have definite values 0 and 1, and the second two are variable, a ∈ {0, 1} and b ∈ {0, 1}, as illustrated in figure 11.These four sites, (0, 1, a, b), influence the values of sites (of different parity) at the next iteration.We consider two cases: the case where the occupied region expands and the sequence at the next iteration is (1, a ′ , b, ), and also the contracting case where the occupations become (0, 1, b ′ ).
In the contracting case, one calculates the probability of b ′ , assuming that a and b are independent, with probability f 1 of being equal to 1: where P 01| (b ′ |a, b) is the conditional probability for obtaining (0, 1, b ′ ) at iteration j + 1 given (0, 1, a, b) at iteration j.Noting that the probability for a 'contracting' shift of the boundary is (1 − p)/2, if the joint distribution of occupations is self-reproducing under iteration, we should expect that Similarly, for the extending case, we can calculate the joint probability distribution of (a ′ , b ′ ): where There does not appear to be any transparent general expression for the conditional probabilities Table 1: Conditional probabilities for reaching states (1, a ′ , b ′ ) or (0, 1, b ′ ) (rows) from initial states (0, (1,0,0) Using the expressions in table 1, we were able to verify that equations (51) and (50) are indeed equal, as are equations ( 53) and (52).We are not aware of any simpler route to this conclusion.
The same argument can be applied at the right-hand edge of the occupied region, so that the site occupations remain independent inside any occupied region, independent of how its boundaries fluctuate.

Relation to standard model for directed percolation
Here we discuss what happens, in one dimension, when p 0 ̸ = 0.This includes the standard model for directed percolation as a special case, and we shall give emphasis to making connections with that problem.So in this section we will not consider the probability distribution of the flux ϕ, but rather concentrate on the set of occupied sites.The model then has just two relevant parameters, p 0 and p 2 .
When p 0 ̸ = 0, there is a possibility for clusters of pathways to die out, so that the occupation fraction f 1 (p 0 , p 2 ) is equal to zero when p 0 > p c , where p c is the percolation threshold.This percolation threshold is a function of p 2 .
The standard directed bond-percolation problem is when the two bonds are filled independently with probability p, so that p 2 ≡ p 2 and (1 − p) 2 = p 0 .The standard directed percolation problem is therefore represented by the parametric line p 0 = (1−p) 2 , p 2 = p 2 in the two-dimensional parameter space of our model.
In sub-section 7.1 we propose a bound on the percolation threshold, p c (p 2 ), and compare this with numerical estimates.Sub-section 7.2 presents some numerical investigations of the critical exponents of the model.We find that these appear to be identical to those of the standard directed percolation model, in accord with a hypothesis of Janssen ( [5]) and Grassberger ([6]).In sub-section 7.3 we present numerical investigations of the PDF of the distribution of sizes of voids.

Phase diagram
The space of models is illustrated in figure 12.The allowed region of p 0 , p 2 space is a triangle, p 0 ≥ 0, p 2 ≥ 0, p 0 + p 2 ≤ 1.The line p 0 = 0 is exactly solvable for the equilibrium distribution as described by equation ( 12) and section 6.The standard directed bond-percolation problem with probability p for bond occupation is the line p 0 = (1 − p) 2 , p 2 = p 2 .We have plotted our numerical evaluation of the critical line, above which f 1 (p 0 , p 2 ) = 0.The critical line crosses the line defining bond-directed percolation at p = 0.644 . .., as expected [12].
We were able to suggest an upper bound on the critical line using the following argument (which assumes that p 2 ≪ 1).We know that p 0 = 0 is exactly solvable, with those sites which are accessible on a given iteration are independently occupied with probability f 1 (0, p 2 ) = 4p 2 /(1 + p 2 ) 2 .Because there are no correlations, the site occupation is a Poisson process.This expression implies that there are 'voids' between occupied sites which have a characteristic lengthscale ⟨L v ⟩ which diverges as p 2 → 0. Assuming that the lattice spacing is unity, and noting that it is only every second site which is accessible, the mean length of the voids is when p 0 = 0.These voids have an exponential distribution of lengths.
A completely empty state is another possible solution.The dynamics of the system is described by the boundary between the empty state and an occupied region.This boundary is described by a single trajectory, and its treatment is more tractable than analysing the joint statistics of an occupied region.The path of the boundary is a random walk with a drift.In the limit where both p 0 and p 2 are small, we can characterise the motion of the boundary of the occupied region by a diffusion coefficient D b and a drift velocity v b .The diffusion coefficient is determined by writing ⟨∆x 2 ⟩ = 2D∆t, and noting that the displacement is unity (except in the rare cases where the trajectory terminates).There is a drift velocity into the empty region, which is determined by noting that when a trajectory splits, with probability p 2 , the daughter trajectory on the void side become the new boundary.When p 0 = 0 the diffusion coefficient D b and a drift velocity v b (into the empty region) are.
Note that after coarse-graining the spatial and temporal scales, the edge of a void satisfies a stochastic differential equation of the form where dη is a standard stochastic increment, satisfying ⟨dη(t)⟩ = 0 and ⟨dη(t)dη(t ′ )⟩ = δ(t − t ′ ).This is equivalent to a Fokker-Planck equation for the position of a void boundary: A similar expression holds for the overall width of the void.Taking account of the fact that the void has two edges, both the drift velocity and the diffusion coefficient are doubled.The steady-state probability density for the distribution of large void sizes is then where the steady state size scale for the large voids is Consider what happens as p 0 is increased.When p 0 ̸ = 0 there is a new mechanism which contributes to the drift velocity of the boundary.When a trail disappears, the boundary of the occupied region retreats by a distance L 1 equal to the length of the first void which is encountered.The velocity of the boundary is then which becomes negative at a critical value of p 0 .The significance of this critical value is that for p 0 > p c the boundary of the occupied region retreats, until we are left with the empty configuration.We should expect that ⟨L 1 ⟩ increases when p 0 > 0, so that (when The critical value of p 0 occurs when v b changes sign, so that In figure 13 we plot p 0 /p 2 2 for the transition line, as a function of p 2 .The ratio is less than 2, as predicted by equation (61).We do not have a theory for the form of this dependence.
Finally, we remark that our model is related to a model of directed percolation in which bonds are deleted with probability 1−p bond , and sites are deleted with probability 1−p site , as discussed in [16].The parameters p bond and p site are related to our parameters as follows: This bond/site deletion model does not cover the whole parameter space of our model: the square 0 ≤ p site ≤ 1, 0 ≤ P bond ≤ 1 maps to the region to the right of the cyan line in figure 12.It can be verified that the region of the critical line of our model lying within this region is in agreement with the data in [16].

Critical exponents
The percolation process which we consider is a form of directed percolation.Grassberger [6] and Janssen [5] introduced a hypothesis that the directed percolation has universal critical exponents, and we expect that the transition in our model lies in this universality class.Thus we expect that f 1 (p 0 , p 2 ) = 0 when p 0 > p c (p 2 ), and that when p c − p 0 is small and positive we have where β = 0.276 . . . is the critical exponent for the order parameter of the directed percolation transition.Figure 14 presents evidence that the occupation fraction vanishes in accord with equation (63).

Void size distribution
We investigated the distribution of sizes of voids, P void (n).The distribution, illustrated in figure 15 is highly distinctive: the is a 'core' region, in which P void (n) has a rapid exponential decay, and a 'tail', which has a slower exponential decay, described by an exponent µ: 2 ), as indicated by the dashed line in the inset.
The lengthscale associated with the slow decay diverges as the critical point is approached.(In figure 15 the void size n is the number of potentially filled sites which are actually empty.Because sites with different parity from the iteration index are automatically empty, the void sizes disccused in sub-section 7.1 are approximately 2n.) Figure 16 shows how values of the decay rate µ behaves as a function of the  16 shows µ as a function of ⟨f ⟩ ν ⊥ /β .Additionally, we divided the value of ⟨f ⟩ ν ⊥ /β by p ν ⊥ c , which reduces the various values to a unique curve.The dashed line shows a power law with a power 1 (linear dependence), which closely agrees with the measured values of µ collapsed to a single curve.

Conclusions, implications for elution
We introduced an alternative model for directed percolation which comes closer to addressing the questions originally posed by Broadbent and Hammersley [1].Our model considers the distribution of fluxes ϕ through wetted bonds, and it has a regime, in which the total flux remains constant.Our model has similarities the Scheidegger river model [18], which also involves addition of fluxes upon combining paths, and which also leads to a power-law distribution of fluxes [19].The models differ as to the source of the flux: in our model the flux total flux is constant when p 0 = 0, whereas in the Scheidegger model new sources are added continuously.
We find that the distribution of fluxes ϕ is very broad, and in section 3 we showed that it has a power-law asymptote at small fluxes: P (ϕ) ∼ ϕ −α .The exponent of this power law, α, was found to vary continuously as a function of the parameters of the model.We found an exact equation for α (equation ( 8)) in terms of some probabilities P k (defined by ( 5)) which characterise the 'wetting' of the skeleton of bonds through which the fluxes flow.
In one dimension, the sites are independently occupied in the steady-state, implying that many quantities of interest, including the exponent α, and the filling fraction for occupied sites, f 1 , can be determined exactly.
In two dimensions, the independent-occupancy assumption gives a very good approximation for α and for the filling fractions, as shown in sections 4 and 5, but it is not exact.In section 6 we examined the conditions for the site occupations to remain independent under iteration.We found that these are only consistent with the predicted values of the filling fraction f in the one-dimensional case.
In the case where the paths can terminate, the model is no longer flux conserving, and we find that the distribution of wetted bonds has a percolation transition.In section 7 we showed that, in the one-dimensional case, the critical exponents are in agreement with those of the standard directed percolation model.We also investigated the distribution of void sizes.
The long-tailed distribution of fluxes is expected to have practical consequences.Consider the following model for an elution process.For definiteness we discuss a solute being washed out of a solid substrate by a liquid.Leaching of a salt from a permeable rock or from land reclaimed from the sea would be a concrete example.We assume that the solid medium is permeated by randomly arranged narrow channels, but that on a large scale it appears homogeneous.We shall assume that the fluid is being forced through the medium, by gravity, or a pressure difference (or both).
If the solid contains solute with a concentration c (defined as mass per unit volume), this will come into equilibrium with the salt dissolved in the liquid phase at a concentration Kc, where K is a partition coefficient.If the liquid is flowing in very narrow channels, we can assume that the solute equilibrates between the solid and liquid phases, so that the concentration in the fluid flowing through the pore is also Kc.The rate at which solute is removed from the medium by flow through the pore is therefore ṁ = Kcϕ, where ϕ is the volume flux through the pore.The concentration is proportional to the amount of solute remaining in the solid phase, so that we expect that the concentration in the runoff through the pore reduces exponentially as a function of time.We therefore expect that the rate of loss of solute from a single pore is where ν is dependent upon the geometry of the pore and the coefficient of solubility in the liquid and solid phases.
In the case of perfusion through a random medium, the liquid may follow many different channels (labelled by an index i), with very different volume fluxes ϕ i in different channels.In this case, the total rate of solute out of the medium is Ṁ (t) = i ṁi (t) = Kc(0) i ϕ i exp(−ν i ϕ i t) . (66) We shall argue that the predominant factor determining the rate of elution at long times is the existence of pores which carry very low fluxes.For this reason we adopt the simplifying assumption that the coefficients ν i are the same for all of the pores.If the probability density function of ϕ is P (ϕ), and the density of pores carrying fluid is ρ, then the time-dependence of the eluted flux from a surface of area A is Ṁ (t) = KρAc(0) where P (s) is the Laplace transform of P (ϕ).Determining the long-time behaviour depends upon the distribution of ϕ in the limit as ϕ → 0. We have argued that this has a power-law behaviour for a wide range of models.This indicates that there are many channels which have an extremely small flux, corresponding to a slow elution of the solute.Using (67), this corresponds a power-law decay of the eluted solute flux: It is noteworthy that, due to the power-law distribution of the flux ϕ, despite the fact that the elution from each channel decreases exponentially as a function of time, the overall rate of elution has a much slower, power-law, decay.This is a consequence of the log-time behaviour being dominated by the channels with the smallest flux.

Figure 4 :
Figure 4: (a) Comparing simulated occupation fraction f 2,2 for the two-branch model with prediction from independent-occupancy approximation, equations (17) and (18).(b) Shows the fractional error, plotted against ln p. (We defined the fractional error by ∆f /f ≡ (f th − f num )/ √ f th f num , where f th and f num are, respectively, the theoretical and numerically determined values of f .)

Figure 5 :
Figure 5: (a) Comparing simulated occupation probability f 2,4 for the four-branch model with prediction from independent-occupancy approximation, equation (21).(b) Shows the fractional error (defined in the same way as for figure4), plotted against ln p.

Figure 6 :
Figure6: Estimates of the probabilities defined in (5), P k (p), k = 1, . . ., 5, for the twodimensional two-branch model.The solid lines are the theoretical predictions, equations (29) to (33).The numerical simulations imposed the requirement that the site has not recently experienced a coalescence in the preceding N thr = 8 iterations.They differ significantly from the theoretical model as p → 1.

results in lim p→0 F
Figure7provides direct tests of the theoretical prediction for the exponent α.We estimated α from simulations with a wide range of values of r and p.In the left panel of Fig.7, we compare the theoretical values of α obtained from equation (9) against numerically estimates, obtained by directly simulating the model: here we used the values of P k obtained from simulations with N thr = 8.In the right panel we collapse all of the data points onto two plots of F (p) ≡ (1 − P 1 )/Q, one obtained using equations (29) to (33), the other using values of P k obtained from simulations with N thr = 8.The former shows small but significant deviations at larger values of p.The small dispersion between the symbols gives an indication of the accuracy of our determination of the exponents α.We remark that, if the errors in the approximations underlying equations (29)-(33) and (17), (18) are negligible as p → 0, we can determine the limiting value of F (p) as p → 0. Using these expressions in equation(10) results in lim p→0 F (p) = 3 .(34)

Figure 7 :
Figure 7: Testing the determination of α for the two-branch two-dimensional, twobranch model, using equation (8), and the P k derived from simulations with N thr = 8.Left panel: there is satisfactory agreement between the empirical values of α and the values obtained from (8).Right panel: the data points collapse onto a plot of F (p) ≡ (1 − P 1 )/Q, using values of P k derived from simulations (with N thr = 8), compared with the values of F (p) obtained from equations (29)-(33).

Figure 8 :
Figure 8: Numerical simulations of the functions P k (p), k = 1, . . ., 7, for the twodimensional, four-branch model.Two different numerical simulations are shown: one includes all sites (i.e.N thr = 0, solid line), the other (points) imposes the requirement that the site has not experienced a coalescence in the preceding N thr = 6 iterations.

Figure 11 :
Figure 11:  We consider configurations at the left-hand edge of the filled region.At iteration j (blue sites), the configuration is (0, 1, a, b) (where a ∈ {0, 1} and b ∈ {0, 1} are variable).At iteration j + 1 (red sites), the occupied region either expands (lower row, occupancy (1, a ′ , b ′ )) or else contracts (upper row, occupancy (0, 1, b ′ )).We show that if a, b are independent, with occupation probability f 1 = 4p/(1 + p) 2 , then a ′ , b ′ have the same property.Sites defined to be empty are shown as open circles, those defined filled are solid colour, and sites which are occupied with probability f 1 (p) are shaded.
P 1| (a, , b ′ |a, b) and P 01| (b ′ |a, b), and we determined them on a case-by-case basis.They are tabulated in table 1.The first four rows determine P 1| (a ′ , b ′ |a, b), and the final two rows specify P 01| (b ′ |a, b).

Figure 12 :
Figure 12: Parameter space of the one-dimensional model, showing the critical line for percolation transition (crosses).The conventional one-dimensional directed percolation problem, in which bonds are occupied independently with probability p, is represented parametrically by the line p 0 = (1 − p) 2 , p 2 = p 2 , shown in blue.The red cross indicates the critical point for the standard directed bond-percolation model, at p ≈ 0.644.

Figure 14 :Figure 15 :
Figure14: Occupied fraction approaching the percolation transition.The average fraction, to the power 1/β, where β is the exponent of directed percolation, shown as a function of p 0 , reveals an approximately linear law.This is consistent with the expected behaviour ⟨f ⟩ ∝ (p c (p 2 ) − p 0 ) β when p 0 → p c (p 2 ), see Eq. (63).

Figure 16 :
Figure16: The decay rate µ of the void size distributions as a function of the normalised mean occupation fraction, ⟨f ⟩/p 2 , to the power ν ⊥ /β, where ν ⊥ and β are the classical directed percolation exponents[3].The dashed line indicate a power 1 corresponds to the directed percolation scaling, and describes approximately the numerical observations. .
2, which is in agreement with equation (54).