# A Size-Perimeter Discrete Growth Model for Percolation Clusters

Authors: Bendegúz Dezső Bak; Tamás Kalmár-Nagy
Journal: Complexity
 (2021)

---

## Abstract

Cluster growth models are utilized for a wide range of scientific and engineering applications, including modeling epidemics and the dynamics of liquid propagation in porous media. Invasion percolation is a stochastic branching process in which a network of sites is getting occupied that leads to the formation of clusters (group of interconnected, occupied sites). The occupation of sites is governed by their resistance distribution; the invasion annexes the sites with the least resistance. An iterative cluster growth model is considered for computing the expected size and perimeter of the growing cluster. A necessary ingredient of the model is the description of the mean perimeter as the function of the cluster size. We propose such a relationship for the site square lattice. The proposed model exhibits (by design) the expected phase transition of percolation models, i.e., it diverges at the percolation thresholdpc. We describe an application for the porosimetry percolation model. The calculations of the cluster growth model compare well with simulation results.


---

## Body

## 1. Introduction

Percolation theory [1–3] has been developed to study the properties of connected clusters in graphs and their associated percolation processes. The simplest problem arising from percolation theory is the site/bond percolation: a regular lattice is considered in which either the cells (sites) or the edges (bonds) are the relevant entities. Each site/bond is independently open with probability p and closed otherwise. The term cluster refers to a set of neighboring open sites/bonds. A fundamental question of percolation theory is whether for a given p, there exists (in the almost sure sense) a cluster that spans through the entire lattice. The probability of the existence of a spanning cluster depends on the occupation probability p and the type of the lattice [4–6]. It was found that there is a strict critical percolation threshold pc characteristic of the lattice type. For an infinite lattice, the probability of the existence of the spanning cluster is zero for p<pc, while it is one for p>pc. At p=pc, there is a singularity in percolation and many properties of the largest cluster follow a power law of p−pc. To this date, pc was derived or computed for many different lattice types. For a few special lattices, the exact value of pc is known, e.g., site percolation threshold of triangular and bond thresholds of triangular, square, and honeycomb lattices [7, 8]. There also exist some other more exotic lattice types, such as the martini lattices and bowtie lattices, for which pc is solved exactly [9]. For other lattice types, an extensive number of simulations were carried out to determine pc numerically [6, 10].The existence of the critical thresholdpc makes percolation suitable to model numerous natural and engineering phenomena [11]. An oft-cited application is the modeling of liquid propagation in a porous medium, while a nowadays highly relevant application is the modeling of disease spread [12, 13]. In the latter, the underlying graph is essentially a network of contacts where the adjacency of the sites defines the contacts. This relates the mean coordination number of the graph, the percolation threshold pc, and the basic reproduction rate (BRP) of diseases that show how many people can be infected by one infected person on average. The critical value of the BRP is 1 and the infection spread becomes unbounded if one person can potentially infect more than one other person. For lattices, the critical BRP is actually larger as each site has only a small group of local contacts; this is why reducing contacts (e.g., travel restrictions) are effective measures of defense. Another way to stop the spread of the disease is to immunize people. Translating this to percolation, the goal is to contain the cluster sizes; i.e., a site is susceptible to infection with probability p and has “immunity” otherwise. Similar to the concept of herd immunity, we can contain the finite cluster sizes, if p<pc [14].We focus on the site variant of the problem, since it is more general than the bond variant (see Section4.2 of [15]). The more complex percolation models were first developed to capture the dynamics of liquid propagation in porous media. The classical invasion percolation model [16, 17] assigns so-called invasion resistances r∈0,1 to the sites of the lattice. The liquid propagation is initiated from a designated set of sites called the starting set. The unoccupied sites that are neighbors of some occupied sites constitute the interface of the occupied cluster (perimeter sites of the cluster). In each step, the liquid occupies the interfacial site whose invasion resistance r is the smallest.Invasion percolation and its variants are complex branching processes [18]. The evolution of the occupied cluster is strongly influenced by its previous states; i.e., the invasion process has memory. Simulation of the invasion process requires tracking of the new interfacial sites that are produced by the recently annexed sites. Our major contribution in this article is the introduction of a simple, iterative cluster growth model to calculate the evolution of the cluster size and perimeter in invasion percolation models. This cluster growth model shares many similarities with the standard epidemic growth algorithm [10, 19].Bak and Kalmár-Nagy introduced the porosimetry percolation model [20, 21], which is a variant of the classical invasion percolation model, to capture the dynamics of intrusion of nonwetting liquid into porous medium, such as in mercury injection porosimetry [22, 23]. In this model, the occupation of the sites is controlled by an external field, the pressurep∈0,1 (analogous to the occupation probability of the classical percolation problem). As opposed to invasion percolation, in porosimetry percolation, multiple sites can be added to the occupied cluster in one step. For a given pressure, every unoccupied, accessible site with r≤p will be occupied. A site is accessible only if it is part of a path of interconnected sites in which every site has r≤p and this path is connected to the starting set.In porosimetry percolation, we can track the evolution of the cluster asp is increased, and compute, e.g., the cluster size as the function of p. In [20], we argued that this percolation simulation is essentially an input-output mapping. The input is the resistance distribution of the sites (cumulative distribution function of the resistance) and the output is cluster size per lattice site as a function of p. Indeed, there is a unique mapping between the input and the output that only depends on the topology of the lattice porosimetry percolation is simulated on. However, a deeper study and explanation were omitted. If this mapping is indeed independent of the input resistance distribution of the sites, is it possible to determine the mapping without simulation? If so, we could use the mapping to predict the cluster size as function of p. Can we also use the cluster growth model to calculate this expected cluster size for any p? In this work, we intend to answer these questions. We show that the considered cluster growth model can predict the cluster size evolution of porosimetry percolation with good accuracy below the percolation threshold pc.This article is structured as follows. In Section2, we summarize the preliminaries: properties of lattice animals and invasion percolation clusters. In Section 3, the iterative cluster growth model is introduced. We derive how the occupied cluster evolves. We also demonstrate that the cluster growth model exhibits the criticality feature of percolation. Details of the cluster growth model are discussed for the site square lattice. In Section 4, we apply the cluster growth model to porosimetry percolation. A detailed description of porosimetry percolation is provided. We prove that the mapping of porosimetry percolation only depends on the topology of the network. We compare the cluster size evolutions obtained with both the cluster growth model and porosimetry percolation simulations. In Section 5, conclusions are drawn.


## 2. Preliminaries

In this section, we define the key terms and the literature results regarding percolation clusters that are the most relevant to this article. The occupied cluster (or simply cluster) is a nearest-neighbor connected setN of occupied sites in ℤ2, also called a lattice animal or polyomino. The size of the polyomino is the cluster size n=N.Another important characteristic is theinterface of the cluster. The interface M is the boundary of N, consisting of the perimeter sites, i.e., the set of sites that are not in N but adjacent to some sites in N. The perimeter is the number of interfacial sites (or perimeter sites), i.e., m=M.Letsn,m denote the number of lattice animals with size n and perimeter m. The so-called perimeter polynomial Dnq is the generating function:(1)Dnq=∑msn,mqm.For example, the first few perimeter polynomials for fixed polyominoes are as follows [24]:(2)D1q=q4,D2q=2q6,D3q=4q7+2q8,D4q=9q8+8q9+2q10.The perimeter polynomials up ton=48 are published on Mertens’ webpage [25]. The number of lattice animals of size n on the site square lattice is(3)sn=∑msn,m=Dn1.The relation between the cluster sizen and the perimeter m is discussed thoroughly in the literature as the statistics of lattice animals [26–28]. The mean perimeter (or mean interface size) of lattice animals of size n is(4)I=1sn∑mmsn,m.As the enumeration ofsn,m is still an open problem [26, 29], I cannot be computed directly from (4). There exist analytic formulae for both the minimum and maximum possible perimeter of a polyomino on the site square lattice as the function of its size [30–32]:(5)minm=2+8n−4,(6)maxm=2n+2.For the mean interface sizeI, there is an empirical formula for large polyominoes proposed by Conway and Guttman [33]:(7)I≈1.195n.Note that lattice animals are purely geometric constructs. In invasion percolation models, the formation of the cluster is governed by the resistances of the sites. A fundamental question is the expected cluster sizeΦ (and interface size I) for a prescribed pressure p such that any site with r>p cannot be occupied.The probability that a cluster of sizen develops from a single starting site is given by the following [34, 35]:(8)∑msn,mpn1−pm.This equation is the link between percolation clusters and lattice animals. Based on (8), one could derive Φ (and I) as the function of p for invasion percolation. However, this derivation again requires the enumeration of lattice animals. There exist a few formulae that describe the expected cluster size Φ as a function of p. Though invasion percolation is different from ordinary percolation, it also exhibits the same criticality feature. This means that above the critical pressure pc, that is equivalent to the percolation threshold, of the invasion becomes unbounded. For p<pc, the occupied cluster consists of a finite number of sites, but for p>pc, it becomes infinitely large. In particular, the percolation threshold for the site square lattice is pc≈0.5927 [36, 37].Researchers were primarily interested in the properties of the cluster close to criticality. Their derivations assumed that the invasion is initiated from a single site of the lattice. Sykes et al. [24] derived the expected cluster size Φ for p’s sufficiently close to pc; that is,(9)Φ≈0.147p−pc−γ,γ=4318.The “low-density” (smallp) series expansion of the expected cluster size Φ as a function of p is in the following form (see, for example, [24]):(10)Φ=∑icipi,i=0,1,…,where ci are constant coefficients. We determined these constants up to c47 based on the perimeter polynomials.In Figure1, we depict the expected cluster size as a function of p for the site square lattice. For the simulations, a finite 1000×1000 lattice was created, and the starting set was an innermost site.Figure 1
Simulated mean cluster size and expected cluster sizes based on Sykes’ formula (9), their low-density series expansion approximation (10) of Φ, and the [46/1] Padé approximant of (10).In this figure, we also depictedΦ based on both (9) and (10). As Figure 1 shows, (9) is only accurate near the criticality, while the series expansion (10) fits well for small p’s but diverges from the simulated curve around p≈0.45. This is why we also calculated the Padé approximant [38] of order 46/1 of (10) because it is accurate on a broad pressure range.There also exists a formula similar to (9) for the mean interface size I of large clusters and p’s sufficiently close to pc that was derived by Stauffer [34] and Wolff and Stauffer [39]. They suggested the following relation:(11)I≈1−ppΦ+1−pp2φxΦσ,σ=3691,where x=p−pcΦσ and φx=−f′x/fx, and(12)fx≈4.81e−7.37x+0.4622.is the scaling function [39]. We can also compute Φ and I from cluster size and interface size statistics of percolation clusters. We carried out 235000 numerical simulations on 1000×1000 site square lattices and extracted the mean interface size I as a function of the mean cluster size Φ (the mean simulation results are analogous to the expected value). The starting set was a single, innermost site in the lattice. The cluster sizes and the interface sizes were evaluated at p=0.5,0.51,…,0.59. Figure 2 shows the mean interface sizes corresponding to different cluster size intervals that are represented by their mean cluster size.Figure 2
The mean interface sizeI of percolation clusters for different cluster size intervals (circles). The thin line connects the first Φ,I=1,4, which is the data point corresponding to a single-site cluster, and the last available data point to highlight the nonlinear relation between I and Φ.
(a)(b)Note that typical large percolation clusters do not obey equation (7); their mean interface size is actually significantly smaller. They rather tend to maintain a “spherical” shape as they spread isotropically on average. Although Figure 2(a) suggests a predominantly linear relation between I and Φ, there is also a softening nonlinear component that is noticeable at smaller Φ’s, see Figure 2(b).Figure3 demonstrates the difference between the perimeter distribution of all possible polyominoes and simulated percolation clusters of the same size.Figure 3
Probability density function of perimeters for all the possible polyominoes of sizen=48 and simulated percolation clusters having the same size.


## 3. The Cluster Growth Model

We consider an iterative cluster growth model to calculate the evolution of percolation clusters. The inputs of the model are the pressurep, the starting cluster size Φ0 and interface size I0 of the prescribed starting set, and the function I=gΦ. We want to calculate the evolution of the expected cluster size Φ and the corresponding mean interface size I. The cluster growth model exhibits the criticality feature of percolation models as the iteration diverges at and above the percolation threshold.We consider uniform resistance distribution of the sites (without the loss of generality, see Section4.1) and a connected starting set. The expected number of sites that get occupied from the initial interface is pI0 at pressure p; this leaves 1−pI0 interfacial sites unoccupied. After the first round of iteration, Φ1 is obtained immediately as(13)Φ1=Φ0+pI0.We havepI0 recently occupied sites; these sites are the “surviving branches.” Each of these sites provides some (e.g., 0–3 for the site square lattice) “new” interfacial sites; let us denote the total number of these new interfacial sites with ΔI. In the next step, the cluster can only occupy the new interfacial sites; thus, the expected number of sites annexed by the cluster is pΔI. Figure 4 explains how the cluster and its interface evolve in one iteration and how to interpret ΔI.Figure 4
(a) The cluster (filled squares) and its interface (dots) after thek−2th iteration step: Φk−2 is the number of black sites and Ik−2 is the number of dots. (b) The cluster (filled and striped squares) and its interface (dots and circles) after the k−1th iteration step: Φk−1 is the number of black and striped sites and Ik−1 is the number of dots and circles. The striped sites are the recently occupied sites, while the circles designate the new interface sites whose number is ΔI. The dots are now “extinct” interfacial sites.
(a)(b)Let us assume thatΦk−1 and Ik−1 are already obtained as in Figure 4(b). The following sites constitute Ik−1:(i)
All the “extinct” interfacial sites that already contributed toIk−2 (dots in Figure 4(b))(ii)
“New” interfacial sites that became part of the interface in thek−1th iteration (circles in Figure 4(b))The extinct interface sites were already tested in thek−1th step. Since they were not occupied, their resistance must be larger than p. Therefore, in the kth step, only the new interface sites will be considered for possible occupation. Based on Figure 4(b), ΔI can be derived and it is(14)ΔI=Ik−1−Ik−2+Φk−1−Φk−2.Thus, the expected number of sites annexed in thekth iteration step is pIk−1−Ik−2+Φk−1−Φk−2. Now, we have Φk and from this train of thought, a simple recursive formula is derived for the expected cluster size and interface size:(15)Φk=Φk−1+pIk−1−Ik−2+Φk−1−Φk−2,k≥2,Ik=gΦk,k≥1,where gΦk is a function of the cluster size. Based on the preliminaries (see Section 2), we expect that gΦk consists of a linear and a nonlinear part, i.e.,(16)Ik=gΦk=glΦk+gnlΦk,glΦk=aΦk+b.We also expectgΦk to be a strictly monotonously increasing function. In Section 3.2, we dissect gΦk further for the site square lattice. Finally, the complete formula of the iteration describing the cluster growth model is as follows:(17)Φ1=Φ0+pI0,Φk=Φk−1+pIk−1−Ik−2+Φk−1−Φk−2,k≥2,Ik=aΦk+b+gnlΦk,k≥1.


### 3.1. Critical Behavior of the Cluster Growth Model

We demonstrate that the cluster growth model exhibits the criticality feature of percolation. The iterative cluster growth calculation yields a difference equation for the cluster size (assuming thatI0=gΦ0 holds); i.e.(18)Φ1=Φ0+pI0,Φk=1+p+paΦk−1−p+paΦk−2+pgnlΦk−1−gnlΦk−2,k≥2.Equation (18) can be transformed into a set of two first-order nonlinear difference equations as follows:(19)Φk−1Φk=01−p+pa1+p+paΦk−2Φk−1+0pgnlΦk−1−gnlΦk−2.Let us investigate the stability of (19) via fixed-point analysis. Linearization of (19) around the fixed point Φ∗ of the iteration yields(20)Φk−1Φk=01−p−pa−pgnl′Φ∗1+p+pa+pgnl′Φ∗Φk−2Φk−1=AΦk−2Φk−1.The eigenvaluesλ1,λ2 of A are obtained by solving the following characteristic equation:(21)λ2−TrAλ+detA=λ2−1+p1+a+gnl′Φ∗λ+p1+a+gnl′Φ∗=0,λ1=1,λ2=p1+a+gnl′Φ∗.Based on Figure2, we expect a softening nonlinearity; i.e., we have(22)limΦ∗⟶∞gnl′Φ∗=0.We get the eigenvaluesλ1=1, λ2=p+pa for Φ∗⟶∞. Since λ1=1, the fixed point is a neutral point for λ2<1 (Lyapunov stable). This latter can only hold for p+pa≤1, so this is a condition for marginal stability. For p+pa>1, the second eigenvalue λ2>1 regardless of the value of Φ∗ and the system becomes unstable. This shows that the radius of convergence of the system (20) is p+pa=1 for large Φ’s. The cluster growth model is consistent with the existence of the critical threshold pc, if(23)a=1−pcpc,holds near pc. This result is also consistent with equation (11).


### 3.2. Interface Update for the Site Square Lattice

We attempt to construct a functionI=gΦ based on literature results that provides accurate calculations compared to percolation simulations. We present here our proposal for the site square lattice, but such a relation can be derived for any other types of lattices.The mean interface sizeI of the cluster of expected size Φ can be calculated with Stauffer’s formula (11). Substituting equation (12) into φx and then φx into (11) leads to(24)I≈1−ppΦ+6.811−pp2Φσ+14.741−pp−1−pp2pcΦ2σ,p<pc.This equation also shows softening nonlinearity as the exponentsσ<2σ<1. We recall that this formula is only accurate for large Φ’s near pc. To make use of (24), we have to determine the lower threshold Φth above which it is applicable with reasonable accuracy. Since we argued that the relation between the mean interface size and the cluster size is softening nonlinear, we expect that the ratio I/Φ decreases for increasing Φ. We can test this expectation by substituting either (9), or the series expansion approximation (10), or the Padé approximant of (10) into (24). We could not work with (9) and (10) due to their inaccuracy in the important midrange p∈0.45,0.55; see Figure 1. However, the Padé approximant is accurate in this particular range; it is a significant improvement over the original series expansion approximation (10). We substituted the Padé approximant into equation (24) and plotted I/Φ against p.Figure5 shows that I/Φ is not decreasing for some p<pth! Via iterative search, we found that I/Φ has a maximum at pth=0.465 (marked by a vertical line in the figure) and the corresponding threshold cluster size is Φth=29.69. This shows that (24) is only applicable for Φ≥Φth=29.69.Figure 5
The ratioI/Φ against p calculated with the Padé approximant of (10). The vertical line marks the maximum of I/Φ on p∈0,pc.Moreover, note that, in the limit ofp⟶pc, equation (24) reduces to(25)I≈1−pcpcΦ+6.811−pcpc2Φσ.We see that the applicability of (24) is heavily restricted: it is only valid for a single-site starting set, Φ0=1 and I0=4. Moreover, we cannot use it, until the expected cluster size exceeds the threshold size Φth=29.69. However, we need to update I also for smaller Φ’s and for different starting sets in the iteration.We assume that below the thresholdΦth, the cluster develops through separate branches. The threshold value Φth=29.69 corresponds to the particular starting set Φ0,ref=1, I0,ref=4 with pth=0.465. We adjust Φth for any larger values of p, Φ0, and I0. A larger Φ0 must increase the threshold size by the same amount, while a larger p and/or I0 increases the number of branches emanating from the starting set by the ratio pI0/pthI0,ref. Equation (26) assigns Φth,mod≥Φth accordingly:(26)Φth,mod=Φ0−Φ0,ref+pI0pthI0,refΦth,Φth=29.69,Φ0,ref=1,I0,ref=4.Let us assume that the starting set is small, i.e.,Φ0≪Φth, but it is not restricted to the single-site case with Φ0=1, I0=4. To calculate Φ and I with high accuracy, a piecewise gΦ is proposed. The linear part glΦ is(27)glΦ=2Φ+2,1≤Φ≤2,0.736Φ+2,2<Φ<Φth,mod,1−ppΦ+2,Φ≥Φth,mod.The nonlinear partgnlΦ is(28)gnlΦ=0,1≤Φ≤2,1.722Φkσ+0.152Φ2σ,2<Φ<Φth,mod,6.811−pp2Φσ+14.741−pp−1−pp2pcΦ2σ,Φ≥Φth,mod.Explanations of the cases are as follows:(i)
For1≤Φ≤2, the expected cluster size Φ is between that of a size of a single site (n=1) and a domino (n=2). In this case, the interface size is trivial; see equation (6): I=maxI=2Φ+2.(ii)
In the range2<Φ<Φth,mod, we use a fit obtained from the mean interface size vs. cluster size statistics from percolation simulations; see Figure 2. The constant parameter b=2 was prescribed to keep consistency with the previous case. The structure of the fit is analogous to Stauffer’s formula but with constant coefficients.(iii)
ForΦ≥Φth,mod, we can finally use Stauffer’s formula (24). The constant parameter b=2 was kept again to keep consistency with the previous case.Calculations with small starting sets are useful to compare the cluster growth model with literature and simulation results. However, for a lot of practical uses, we need very large starting sets for whichΦ0≫Φth. For such large starting sets, we will simply use(29)I=gΦ=1−pcpcΦ+2+6.811−pcpc2Φσ,Φ0≫Φth.That is consistent with the limit expression (25) of Stauffer’s formula for large cluster sizes. The relaxation of Ik is required at some iteration steps to avoid two possible contradictions. For starting sets with I0≫gΦ0, it is possible to get Ik≪Ik−1 for low p’s. This is highly atypical during the low-pressure build-up of percolation clusters. Moreover, we experienced that Stauffer’s formula may yield Ik/Φk>Ik−1/Φk−1 when p≪pc, whereas I/Φ ratio must decrease as Φ increases; see Figure 5. Therefore, we also incorporate a simple relaxation:(30)Ik=Ik−1,gΦk<Ik−1,gΦk,Ik−1≤gΦk≤Ik−1Φk−1Φk,Ik−1Φk−1Φk,gΦk>Ik−1Φk−1Φk.Now, we show an example of howΦ and I evolve in the cluster growth model. We study the Φ,I plane depicted in Figure 6. All the possible Φ,I pairs are located within the white area, which is bounded by the theoretical limits based on equations (5) and (6). In Figure 6, the dots show the evolution of Φk,Ik for pressure pth≤p<pc and starting Φ0≪Φth. By design, the iteration attracts any pair of initial values Φ0,I0 onto I=gΦ in one iteration step. Thus, the iteration jumps to Φ1,I1 located on gΦ and then follows gΦ towards the fixed point Φ∗,I∗. Upon reaching Φth,mod, the function gΦ changes as specified by (27) and (28). Hence, there is a slight but still noticeable breakpoint here. Further breakpoints can be caused by the relaxation (30), which was the case in the depicted example.Figure 6
The evolution of the iteration shown in theΦ,I plane. The initial position and the fixed point are denoted by empty and filled squares, respectively. The solid line and the dots show the evolution of Φk,Ik. The vertical line marks Φth,mod.


## 3.1. Critical Behavior of the Cluster Growth Model

We demonstrate that the cluster growth model exhibits the criticality feature of percolation. The iterative cluster growth calculation yields a difference equation for the cluster size (assuming thatI0=gΦ0 holds); i.e.(18)Φ1=Φ0+pI0,Φk=1+p+paΦk−1−p+paΦk−2+pgnlΦk−1−gnlΦk−2,k≥2.Equation (18) can be transformed into a set of two first-order nonlinear difference equations as follows:(19)Φk−1Φk=01−p+pa1+p+paΦk−2Φk−1+0pgnlΦk−1−gnlΦk−2.Let us investigate the stability of (19) via fixed-point analysis. Linearization of (19) around the fixed point Φ∗ of the iteration yields(20)Φk−1Φk=01−p−pa−pgnl′Φ∗1+p+pa+pgnl′Φ∗Φk−2Φk−1=AΦk−2Φk−1.The eigenvaluesλ1,λ2 of A are obtained by solving the following characteristic equation:(21)λ2−TrAλ+detA=λ2−1+p1+a+gnl′Φ∗λ+p1+a+gnl′Φ∗=0,λ1=1,λ2=p1+a+gnl′Φ∗.Based on Figure2, we expect a softening nonlinearity; i.e., we have(22)limΦ∗⟶∞gnl′Φ∗=0.We get the eigenvaluesλ1=1, λ2=p+pa for Φ∗⟶∞. Since λ1=1, the fixed point is a neutral point for λ2<1 (Lyapunov stable). This latter can only hold for p+pa≤1, so this is a condition for marginal stability. For p+pa>1, the second eigenvalue λ2>1 regardless of the value of Φ∗ and the system becomes unstable. This shows that the radius of convergence of the system (20) is p+pa=1 for large Φ’s. The cluster growth model is consistent with the existence of the critical threshold pc, if(23)a=1−pcpc,holds near pc. This result is also consistent with equation (11).


## 3.2. Interface Update for the Site Square Lattice

We attempt to construct a functionI=gΦ based on literature results that provides accurate calculations compared to percolation simulations. We present here our proposal for the site square lattice, but such a relation can be derived for any other types of lattices.The mean interface sizeI of the cluster of expected size Φ can be calculated with Stauffer’s formula (11). Substituting equation (12) into φx and then φx into (11) leads to(24)I≈1−ppΦ+6.811−pp2Φσ+14.741−pp−1−pp2pcΦ2σ,p<pc.This equation also shows softening nonlinearity as the exponentsσ<2σ<1. We recall that this formula is only accurate for large Φ’s near pc. To make use of (24), we have to determine the lower threshold Φth above which it is applicable with reasonable accuracy. Since we argued that the relation between the mean interface size and the cluster size is softening nonlinear, we expect that the ratio I/Φ decreases for increasing Φ. We can test this expectation by substituting either (9), or the series expansion approximation (10), or the Padé approximant of (10) into (24). We could not work with (9) and (10) due to their inaccuracy in the important midrange p∈0.45,0.55; see Figure 1. However, the Padé approximant is accurate in this particular range; it is a significant improvement over the original series expansion approximation (10). We substituted the Padé approximant into equation (24) and plotted I/Φ against p.Figure5 shows that I/Φ is not decreasing for some p<pth! Via iterative search, we found that I/Φ has a maximum at pth=0.465 (marked by a vertical line in the figure) and the corresponding threshold cluster size is Φth=29.69. This shows that (24) is only applicable for Φ≥Φth=29.69.Figure 5
The ratioI/Φ against p calculated with the Padé approximant of (10). The vertical line marks the maximum of I/Φ on p∈0,pc.Moreover, note that, in the limit ofp⟶pc, equation (24) reduces to(25)I≈1−pcpcΦ+6.811−pcpc2Φσ.We see that the applicability of (24) is heavily restricted: it is only valid for a single-site starting set, Φ0=1 and I0=4. Moreover, we cannot use it, until the expected cluster size exceeds the threshold size Φth=29.69. However, we need to update I also for smaller Φ’s and for different starting sets in the iteration.We assume that below the thresholdΦth, the cluster develops through separate branches. The threshold value Φth=29.69 corresponds to the particular starting set Φ0,ref=1, I0,ref=4 with pth=0.465. We adjust Φth for any larger values of p, Φ0, and I0. A larger Φ0 must increase the threshold size by the same amount, while a larger p and/or I0 increases the number of branches emanating from the starting set by the ratio pI0/pthI0,ref. Equation (26) assigns Φth,mod≥Φth accordingly:(26)Φth,mod=Φ0−Φ0,ref+pI0pthI0,refΦth,Φth=29.69,Φ0,ref=1,I0,ref=4.Let us assume that the starting set is small, i.e.,Φ0≪Φth, but it is not restricted to the single-site case with Φ0=1, I0=4. To calculate Φ and I with high accuracy, a piecewise gΦ is proposed. The linear part glΦ is(27)glΦ=2Φ+2,1≤Φ≤2,0.736Φ+2,2<Φ<Φth,mod,1−ppΦ+2,Φ≥Φth,mod.The nonlinear partgnlΦ is(28)gnlΦ=0,1≤Φ≤2,1.722Φkσ+0.152Φ2σ,2<Φ<Φth,mod,6.811−pp2Φσ+14.741−pp−1−pp2pcΦ2σ,Φ≥Φth,mod.Explanations of the cases are as follows:(i)
For1≤Φ≤2, the expected cluster size Φ is between that of a size of a single site (n=1) and a domino (n=2). In this case, the interface size is trivial; see equation (6): I=maxI=2Φ+2.(ii)
In the range2<Φ<Φth,mod, we use a fit obtained from the mean interface size vs. cluster size statistics from percolation simulations; see Figure 2. The constant parameter b=2 was prescribed to keep consistency with the previous case. The structure of the fit is analogous to Stauffer’s formula but with constant coefficients.(iii)
ForΦ≥Φth,mod, we can finally use Stauffer’s formula (24). The constant parameter b=2 was kept again to keep consistency with the previous case.Calculations with small starting sets are useful to compare the cluster growth model with literature and simulation results. However, for a lot of practical uses, we need very large starting sets for whichΦ0≫Φth. For such large starting sets, we will simply use(29)I=gΦ=1−pcpcΦ+2+6.811−pcpc2Φσ,Φ0≫Φth.That is consistent with the limit expression (25) of Stauffer’s formula for large cluster sizes. The relaxation of Ik is required at some iteration steps to avoid two possible contradictions. For starting sets with I0≫gΦ0, it is possible to get Ik≪Ik−1 for low p’s. This is highly atypical during the low-pressure build-up of percolation clusters. Moreover, we experienced that Stauffer’s formula may yield Ik/Φk>Ik−1/Φk−1 when p≪pc, whereas I/Φ ratio must decrease as Φ increases; see Figure 5. Therefore, we also incorporate a simple relaxation:(30)Ik=Ik−1,gΦk<Ik−1,gΦk,Ik−1≤gΦk≤Ik−1Φk−1Φk,Ik−1Φk−1Φk,gΦk>Ik−1Φk−1Φk.Now, we show an example of howΦ and I evolve in the cluster growth model. We study the Φ,I plane depicted in Figure 6. All the possible Φ,I pairs are located within the white area, which is bounded by the theoretical limits based on equations (5) and (6). In Figure 6, the dots show the evolution of Φk,Ik for pressure pth≤p<pc and starting Φ0≪Φth. By design, the iteration attracts any pair of initial values Φ0,I0 onto I=gΦ in one iteration step. Thus, the iteration jumps to Φ1,I1 located on gΦ and then follows gΦ towards the fixed point Φ∗,I∗. Upon reaching Φth,mod, the function gΦ changes as specified by (27) and (28). Hence, there is a slight but still noticeable breakpoint here. Further breakpoints can be caused by the relaxation (30), which was the case in the depicted example.Figure 6
The evolution of the iteration shown in theΦ,I plane. The initial position and the fixed point are denoted by empty and filled squares, respectively. The solid line and the dots show the evolution of Φk,Ik. The vertical line marks Φth,mod.


## 4. Application for Porosimetry Percolation

A porosimetry percolation simulation can compute the evolution of the cluster for anyp∈0,1, Φ can be computed for any p, and it is finite on a finite network of sites. The cluster growth model can only calculate the evolution of the cluster for pressures below pc, since the iteration diverges for p≥pc. Therefore, we cannot replace porosimetry percolation simulation with the current cluster growth model. Still, such a model has high practical importance since the larger pores and most of the volume of a porous medium are invaded typically at low pressures. Therefore, geologists and reservoir engineers particularly focus on the low pressure part of measurement results.We note that in the approaches cited below, percolation is simulated mostly on simple, connected graphs that represent three-dimensional pore structure rather than on regular lattices. Site lattices are special subsets of simple, connected graphs, as any site lattice can be defined as a simple, connected graph. The sites are the vertices of the graph, and their adjacency defines the edges of the graph.Our main goal with porosimetry percolation was to exploit the most important result of a mercury injection porosimetry measurement, the saturation curve. The saturationϕ shows the fraction of the total pore volume occupied by mercury at pressure p, the saturation curve shows ϕ as a function of p. In the terminology of percolation and cluster growth models, the saturation can be interpreted as the expected cluster size scaled by the size of the network, e.g.,(31)ϕ=ΦL2,where L is the side length of a finite site square lattice. The Washburn equation [40] specifies the relation p∼1/ρ between pressure p and pore size ρ of the smallest pore, which can be filled up with mercury at that pressure. That is, a larger p is required to fill a pore characterized by a smaller ρ. The pore size ρ is a physical property that is interpreted in percolation models as r∼1/ρ, the resistance of the pore. We clarify that resistance r is a pressure-type quantity; it is equal to pressure p at which the pore can be occupied. Based on these relations, the resistance distribution of pores and thus the pore size distribution (PSD) can be derived from the saturation curve [41].However, this derived PSD is incorrect due to the shielding mechanism [42–44]. Shielding occurs when narrow capillaries surround a large pore, preventing it from being filled up with mercury at the pressure equal to its resistance; hence, this large pore is shielded. It will be occupied at a higher pressure, so it will contribute to the saturation curve at that higher pressure. Therefore, any resistance distribution and PSD derived from the saturation curve are distorted compared to the real PSD.Many studies discuss how to fix the PSD derived from the saturation curve. Iterative optimization methods [44–46] that take the PSD derived from measurement data as an initial guess are a popular approach. The optimization goal is to find a PSD for which the simulated saturation is close to the measured one. There are other methods as well that do not require material injection into the porous sample; e.g., imaging techniques [47–49] are common.Bak and Kalmár-Nagy [20] proposed another method to compute the PSD of porous rock samples, which involves porosimetry percolation simulations on a random graph model of the pore network. Compared to previously mentioned iterative methods, which were similar in spirit, their method does not require an iterative search to get the correct PSD. This method also requires an initial input distribution that is derived from the measured saturation curve that we denote with ϕ1p. As an initial guess, the PSD derived from ϕ1p is used to assign a preliminary resistance distribution to the network on which porosimetry percolation is simulated. The result of the simulation is another saturation curve, denoted by ϕ2p. The main idea was to treat the simulation as an input-output mapping: from input ϕ1p, it yielded output ϕ2p. The task is to find input ϕ0p for which the output is the measured ϕ1p. Then, the resistance distribution and the PSD can be derived from ϕ0p instead of ϕ1p. The method is illustrated in Figure 7.Figure 7
The unknownϕ0p, the experimental saturation curve ϕ1p, and the simulated saturation curve ϕ2p.The mappingP relates ϕ1p and ϕ2p as(32)ϕ2p=Pϕ1p.Since the saturation curves are strictly monotonously increasing functions onp∈0,1, P is a continuous bijection, whose inverse is also continuous; i.e., P is a homeomorphism. Thus, the inverse of P exists and it can be used to obtain ϕ0p as(33)ϕ0p=P−1ϕ1p.It was argued in [20] that P does not depend on the input ϕ1p, it only depends on the structure of the network porosimetry percolation is simulated on. We provide a proof to this claim in the next subsection.


### 4.1. Porosimetry Percolation Mapping

We consider a lattice onℤ2 and the starting set is a single site. Since resistance r is a pressure-type quantity, we can describe the input resistance distribution with a cumulative distribution function (CDF) Fp≡Fr.DistributionFp shows the probability that the resistance of a site is r≤p. Based on equation (8), the probability that the occupied cluster will have a size of n is given by(34)∑msn,mFnp1−Fpm.The expected cluster size is(35)Φ=∑nn∑msn,mFnp1−Fpm.Assuming a finite network of sites, the saturation isΦ scaled with the number of sites in the lattice. For instance, the expected saturation of a site square lattice with L2 sites is as follows:(36)ϕ=1L2∑nn∑msn,mFnp1−Fpm.We can calculate (36) for all p∈0,1 to get the saturation curve ϕp.Let us draw the resistances of the sites from a uniform distribution; i.e., we haveFp=p. In this case, the identity function is transformed into ϕp by the mapping P; i.e.,(37)ϕp=Pp.Do we get the sameϕ for any resistance distribution Fp? Based on (36) and (37), we can formulate(38)ϕp=Pp=1L2∑nn∑msn,mpn1−pm,p∈0,1.SubstitutingFp for p into (38) yields(39)ϕp=PFp=1L2∑nn∑msn,mFnp1−Fpm,p∈0,1.We got back equation (32) on the left-hand side with different notations. This proves that mapping P is indeed independent of Fp.The independence ofP from Fp is an important result since it means we can predict the expected cluster size/saturation as the function of the pressure without actual simulations for any Fp, e.g., using the proposed cluster growth model.


### 4.2. Tests with Small Starting Sets (Φ0≪Φth)

We test the cluster growth model on the site square lattice with the interface update proposed in Section3.2. Based on Section 4.1, we can draw the resistances of the sites from uniform distribution without loss of generality.First, we calculated the evolution ofΦ and I for p∈0,pc for three different starting sets: single site (Φ0=1,I0=4), maximum perimeter tromino (Φ0=3,I0=8), and maximum perimeter heptomino (Φ0=7,I0=16). We also carried out porosimetry percolation simulations with these starting sets on a site square lattice of size L=1000 (the number of sites was L2=106). In the simulations, the sites constituting the starting set were always located in the innermost position within the lattice, as far as possible from its boundaries. In Figure 8, the evolution of the expected cluster size Φk vs. the iteration number k is depicted for the single-site starting set; the parameter of the curves is p.Figure 8
(a) The evolution ofΦk for single-site starting set with different p’s. (b) The evolution of Φk for single-site starting set with p=0.25 and p=0.55. The thick lines are calculated with the cluster growth model and the thin lines are individual simulation results.
(a)(b)Figure8 shows that the evolution calculated with the cluster growth model is close to the simulated evolution for p<pc. Results are also shown for some p≥pc cases to demonstrate that the iteration indeed diverges for those pressures.In Figure9, Φ/Φ0 is shown for the pressure range p∈0,pc with the three different starting sets. In Figure 10, the results for the single-site starting set are also shown together with the expected cluster size calculated with (9). The calculations of the cluster growth model were performed with k=10000 iterations, i.e., Φ=Φ10000. The simulation results are averaged over 1000 simulations.Figure 9
Φ scaled with the starting size Φ0 calculated using both the cluster growth model and porosimetry percolation simulations. The starting sets were (a) single site (Φ0=1, I0=4), (b) maximum perimeter tromino (Φ0=3, I0=8), and (c) maximum perimeter heptomino (Φ0=7, I0=16).
(a)(b)(c)Figure 10
Φ scaled with the starting size Φ0 vs. pc−p on a log-log plot. The starting set was a single site (Φ0=1, I0=4).Figures9 and 10 show that the cluster growth model provides a good approximation of Φ for 0≤p<pc, though the calculation overshoots both the simulation result and Sykes’ expected cluster size close to pc. The calculation is also convincing for the tromino and heptomino starting sets. This justifies the proposed interface update, which is described by equations (27) and (28).To quantify the accuracy of the iterative calculation compared to porosimetry percolation simulation results, the relative error of the iterative cluster growth calculation against the simulation data was computed forp=0.01i, i=0,1,…,j, where j=55,58,59 were chosen. The mean and root mean squared (RMS) of these relative errors are listed in Table 1.Table 1
Mean and RMS of the relative errorsΦcalc−Φsim/Φsim of the cluster growth calculation for small starting sets.
CaseΦcalc−Φsim/Φsim¯ (%)1/j∑Φcalc−Φsim/Φsim2 (%)Φ0=1, I0=4j=555.27.4j=587.614.2j=5917.578.8Φ0=3, I0=8j=556.78.4j=588.010.9j=5916.164.8Φ0=7, I0=16j=558.810.4j=589.010.9j=5916.156.9According to the relative errors, the cluster growth model considerably overshoots the simulation result nearpc. Still, the relative errors are acceptable if we exclude the data points that are very close to pc. In fact, in the simulations, the finite size of the lattice limits the expansion of the cluster. We argue that the simulation result would be closer to the cluster growth model near threshold pc, if the simulations were carried out on a much larger lattice.In Figure11, representative visualizations of the evolution on the Φ,I plane are shown for p=0.55. The evolution of the cluster growth model shows nice agreement with the simulated evolutions.Figure 11
Simulated and calculated evolutions on theΦ,I plane for (a) single site (Φ0=1, I0=4), (b) maximum perimeter tromino (Φ0=3, I0=8), and (c) maximum perimeter heptomino (Φ0=7, I0=16) starting sets. The thick lines show the evolution of the cluster growth model and the thin lines are evolutions computed with individual porosimetry percolation simulations. The empty squares and full squares correspond to the starting and fully evolved clusters, respectively.
(a)(b)(c)We can conclude based on the results that the cluster growth model predicts the evolution of the occupied cluster with good accuracy in the low-pressure range0≤p<pc.


### 4.3. Results for Large Starting Sets (Φ0≫Φth)

Calculation with large initial cluster sizes is of high practical importance. Consider a simulation in which the starting set is all the sites that are located on one or multiple edge boundaries of the lattice. In other words, we choose one or more sides of the lattice to be the starting set. The motivation behind this choice of the starting set is that it can be a simple model for the invasion of a dry rock sample immersed in nonwetting liquid.We use the cluster growth model to predict the evolution ofΦ and I for inputs Φ0,I0 that correspond to the size and perimeter of one or more sides of a finite lattice. That is, for a 1-sided invasion of the site square lattice of size L, we have(40)Φ0=L,I0=L.In the 4-sided invasion case, the starting set completely surrounds the rest of the lattice; we have(41)Φ0=4L−4,I0=4L−2−4.We carried out the calculations with the cluster growth model and porosimetry percolation simulations. We investigated 6 cases altogether: 2 different starting set configurations (1-sided and 4-sided invasions) and 3 different lattice sizes (L=500,1000,2000). In Figure 12, the expected saturation ϕ=Φ/L2 of the lattice is shown in the pressure range 0≤p<pc for the different starting sets. The calculations of the cluster growth model were performed with k=10000 iterations, i.e., Φ=Φ10000. The simulation results are averages of 1000, 500, and 100 simulations for L=500,1000,2000, respectively.Figure 12
ϕ=Φ/L2 (expected cluster size scaled with the lattice size) calculated using both the cluster growth model and porosimetry percolation simulations. 1-sided invasion of lattices having (a) L=500 (Φ0=500, I0=500), (b) L=1000 (Φ0=1000, I0=1000), and (c) L=2000 (Φ0=2000, I0=2000). 4-sided invasion of lattices having (d) L=500 (Φ0=1996, I0=1988), (e) L=1000 (Φ0=3996, I0=3988), and (f) L=2000 (Φ0=7996, I0=7988).
(a)(b)(c)(d)(e)(f)The plots are convincing regarding the accuracy of the cluster growth model. Similar to Section4.2, the relative error of the iterative cluster growth calculation against the simulation data was computed for p=0.01i, i=0,…,j, where j=55,58,59 were chosen. The mean and RMS of these relative errors are listed in Table 2.Table 2
Mean and RMS of the relative errorsΦcalc−Φsim/Φsim of the cluster growth calculation for 1-sided and 4-sided invasion cases.
CaseΦcalc−Φsim/Φsim¯ (%)1/j∑Φcalc−Φsim/Φsim2 (%)1-sidedL=500j=551.42.3j=582.14.3j=596.131.51-sidedL=1000j=552.63.4j=582.83.7j=593.88.51-sidedL=2000j=554.96.8j=585.78.1j=595.78.14-sidedL=500j=554.35.7j=584.56.0j=598.733.24-sidedL=1000j=555.88.1j=586.79.4j=597.813.24-sidedL=2000j=556.89.6j=588.011.4j=597.911.4The listed relative errors demonstrate that the estimation of the cluster size with the cluster growth model is qualitatively good and also quantitatively acceptable for these starting sets. For the smallestL=500 case, the expected cluster size significantly overshoots the mean simulation result near pc. Again, the finite size of the lattice explains the issue: in the simulations, the expansion of the cluster is confined by the lattice boundaries, while there is no such restriction in the cluster growth model. For larger lattices, this becomes a lesser issue.In Figure13, representative visualizations of the evolution on the Φ,I plane are shown for p=0.55.Figure 13
Calculated and simulated evolutions on theΦ,I plane for the following cases: 1-sided invasion of lattices having (a) L=1000 (Φ0=1000, I0=1000), (b) L=2000 (Φ0=2000, I0=2000), and 4-sided invasion of lattices having (c) L=500 (Φ0=1996, I0=1988), (d) L=1000 (Φ0=3996, I0=3988). The thick lines show the evolution of the cluster growth model; the thin lines are evolutions computed with individual porosimetry percolation simulations. The empty squares and full squares correspond to the starting and fully evolved clusters, respectively.
(a)(b)(c)(d)With everything considered, the results prove the viability of the cluster growth model for possible application on predicting the prime characteristics (Φ and I) of the occupied cluster in lattices (or in simple connected graphs) for the low-pressure range 0≤p<pc.


## 4.1. Porosimetry Percolation Mapping

We consider a lattice onℤ2 and the starting set is a single site. Since resistance r is a pressure-type quantity, we can describe the input resistance distribution with a cumulative distribution function (CDF) Fp≡Fr.DistributionFp shows the probability that the resistance of a site is r≤p. Based on equation (8), the probability that the occupied cluster will have a size of n is given by(34)∑msn,mFnp1−Fpm.The expected cluster size is(35)Φ=∑nn∑msn,mFnp1−Fpm.Assuming a finite network of sites, the saturation isΦ scaled with the number of sites in the lattice. For instance, the expected saturation of a site square lattice with L2 sites is as follows:(36)ϕ=1L2∑nn∑msn,mFnp1−Fpm.We can calculate (36) for all p∈0,1 to get the saturation curve ϕp.Let us draw the resistances of the sites from a uniform distribution; i.e., we haveFp=p. In this case, the identity function is transformed into ϕp by the mapping P; i.e.,(37)ϕp=Pp.Do we get the sameϕ for any resistance distribution Fp? Based on (36) and (37), we can formulate(38)ϕp=Pp=1L2∑nn∑msn,mpn1−pm,p∈0,1.SubstitutingFp for p into (38) yields(39)ϕp=PFp=1L2∑nn∑msn,mFnp1−Fpm,p∈0,1.We got back equation (32) on the left-hand side with different notations. This proves that mapping P is indeed independent of Fp.The independence ofP from Fp is an important result since it means we can predict the expected cluster size/saturation as the function of the pressure without actual simulations for any Fp, e.g., using the proposed cluster growth model.


## 4.2. Tests with Small Starting Sets (Φ0≪Φth)

We test the cluster growth model on the site square lattice with the interface update proposed in Section3.2. Based on Section 4.1, we can draw the resistances of the sites from uniform distribution without loss of generality.First, we calculated the evolution ofΦ and I for p∈0,pc for three different starting sets: single site (Φ0=1,I0=4), maximum perimeter tromino (Φ0=3,I0=8), and maximum perimeter heptomino (Φ0=7,I0=16). We also carried out porosimetry percolation simulations with these starting sets on a site square lattice of size L=1000 (the number of sites was L2=106). In the simulations, the sites constituting the starting set were always located in the innermost position within the lattice, as far as possible from its boundaries. In Figure 8, the evolution of the expected cluster size Φk vs. the iteration number k is depicted for the single-site starting set; the parameter of the curves is p.Figure 8
(a) The evolution ofΦk for single-site starting set with different p’s. (b) The evolution of Φk for single-site starting set with p=0.25 and p=0.55. The thick lines are calculated with the cluster growth model and the thin lines are individual simulation results.
(a)(b)Figure8 shows that the evolution calculated with the cluster growth model is close to the simulated evolution for p<pc. Results are also shown for some p≥pc cases to demonstrate that the iteration indeed diverges for those pressures.In Figure9, Φ/Φ0 is shown for the pressure range p∈0,pc with the three different starting sets. In Figure 10, the results for the single-site starting set are also shown together with the expected cluster size calculated with (9). The calculations of the cluster growth model were performed with k=10000 iterations, i.e., Φ=Φ10000. The simulation results are averaged over 1000 simulations.Figure 9
Φ scaled with the starting size Φ0 calculated using both the cluster growth model and porosimetry percolation simulations. The starting sets were (a) single site (Φ0=1, I0=4), (b) maximum perimeter tromino (Φ0=3, I0=8), and (c) maximum perimeter heptomino (Φ0=7, I0=16).
(a)(b)(c)Figure 10
Φ scaled with the starting size Φ0 vs. pc−p on a log-log plot. The starting set was a single site (Φ0=1, I0=4).Figures9 and 10 show that the cluster growth model provides a good approximation of Φ for 0≤p<pc, though the calculation overshoots both the simulation result and Sykes’ expected cluster size close to pc. The calculation is also convincing for the tromino and heptomino starting sets. This justifies the proposed interface update, which is described by equations (27) and (28).To quantify the accuracy of the iterative calculation compared to porosimetry percolation simulation results, the relative error of the iterative cluster growth calculation against the simulation data was computed forp=0.01i, i=0,1,…,j, where j=55,58,59 were chosen. The mean and root mean squared (RMS) of these relative errors are listed in Table 1.Table 1
Mean and RMS of the relative errorsΦcalc−Φsim/Φsim of the cluster growth calculation for small starting sets.
CaseΦcalc−Φsim/Φsim¯ (%)1/j∑Φcalc−Φsim/Φsim2 (%)Φ0=1, I0=4j=555.27.4j=587.614.2j=5917.578.8Φ0=3, I0=8j=556.78.4j=588.010.9j=5916.164.8Φ0=7, I0=16j=558.810.4j=589.010.9j=5916.156.9According to the relative errors, the cluster growth model considerably overshoots the simulation result nearpc. Still, the relative errors are acceptable if we exclude the data points that are very close to pc. In fact, in the simulations, the finite size of the lattice limits the expansion of the cluster. We argue that the simulation result would be closer to the cluster growth model near threshold pc, if the simulations were carried out on a much larger lattice.In Figure11, representative visualizations of the evolution on the Φ,I plane are shown for p=0.55. The evolution of the cluster growth model shows nice agreement with the simulated evolutions.Figure 11
Simulated and calculated evolutions on theΦ,I plane for (a) single site (Φ0=1, I0=4), (b) maximum perimeter tromino (Φ0=3, I0=8), and (c) maximum perimeter heptomino (Φ0=7, I0=16) starting sets. The thick lines show the evolution of the cluster growth model and the thin lines are evolutions computed with individual porosimetry percolation simulations. The empty squares and full squares correspond to the starting and fully evolved clusters, respectively.
(a)(b)(c)We can conclude based on the results that the cluster growth model predicts the evolution of the occupied cluster with good accuracy in the low-pressure range0≤p<pc.


## 4.3. Results for Large Starting Sets (Φ0≫Φth)

Calculation with large initial cluster sizes is of high practical importance. Consider a simulation in which the starting set is all the sites that are located on one or multiple edge boundaries of the lattice. In other words, we choose one or more sides of the lattice to be the starting set. The motivation behind this choice of the starting set is that it can be a simple model for the invasion of a dry rock sample immersed in nonwetting liquid.We use the cluster growth model to predict the evolution ofΦ and I for inputs Φ0,I0 that correspond to the size and perimeter of one or more sides of a finite lattice. That is, for a 1-sided invasion of the site square lattice of size L, we have(40)Φ0=L,I0=L.In the 4-sided invasion case, the starting set completely surrounds the rest of the lattice; we have(41)Φ0=4L−4,I0=4L−2−4.We carried out the calculations with the cluster growth model and porosimetry percolation simulations. We investigated 6 cases altogether: 2 different starting set configurations (1-sided and 4-sided invasions) and 3 different lattice sizes (L=500,1000,2000). In Figure 12, the expected saturation ϕ=Φ/L2 of the lattice is shown in the pressure range 0≤p<pc for the different starting sets. The calculations of the cluster growth model were performed with k=10000 iterations, i.e., Φ=Φ10000. The simulation results are averages of 1000, 500, and 100 simulations for L=500,1000,2000, respectively.Figure 12
ϕ=Φ/L2 (expected cluster size scaled with the lattice size) calculated using both the cluster growth model and porosimetry percolation simulations. 1-sided invasion of lattices having (a) L=500 (Φ0=500, I0=500), (b) L=1000 (Φ0=1000, I0=1000), and (c) L=2000 (Φ0=2000, I0=2000). 4-sided invasion of lattices having (d) L=500 (Φ0=1996, I0=1988), (e) L=1000 (Φ0=3996, I0=3988), and (f) L=2000 (Φ0=7996, I0=7988).
(a)(b)(c)(d)(e)(f)The plots are convincing regarding the accuracy of the cluster growth model. Similar to Section4.2, the relative error of the iterative cluster growth calculation against the simulation data was computed for p=0.01i, i=0,…,j, where j=55,58,59 were chosen. The mean and RMS of these relative errors are listed in Table 2.Table 2
Mean and RMS of the relative errorsΦcalc−Φsim/Φsim of the cluster growth calculation for 1-sided and 4-sided invasion cases.
CaseΦcalc−Φsim/Φsim¯ (%)1/j∑Φcalc−Φsim/Φsim2 (%)1-sidedL=500j=551.42.3j=582.14.3j=596.131.51-sidedL=1000j=552.63.4j=582.83.7j=593.88.51-sidedL=2000j=554.96.8j=585.78.1j=595.78.14-sidedL=500j=554.35.7j=584.56.0j=598.733.24-sidedL=1000j=555.88.1j=586.79.4j=597.813.24-sidedL=2000j=556.89.6j=588.011.4j=597.911.4The listed relative errors demonstrate that the estimation of the cluster size with the cluster growth model is qualitatively good and also quantitatively acceptable for these starting sets. For the smallestL=500 case, the expected cluster size significantly overshoots the mean simulation result near pc. Again, the finite size of the lattice explains the issue: in the simulations, the expansion of the cluster is confined by the lattice boundaries, while there is no such restriction in the cluster growth model. For larger lattices, this becomes a lesser issue.In Figure13, representative visualizations of the evolution on the Φ,I plane are shown for p=0.55.Figure 13
Calculated and simulated evolutions on theΦ,I plane for the following cases: 1-sided invasion of lattices having (a) L=1000 (Φ0=1000, I0=1000), (b) L=2000 (Φ0=2000, I0=2000), and 4-sided invasion of lattices having (c) L=500 (Φ0=1996, I0=1988), (d) L=1000 (Φ0=3996, I0=3988). The thick lines show the evolution of the cluster growth model; the thin lines are evolutions computed with individual porosimetry percolation simulations. The empty squares and full squares correspond to the starting and fully evolved clusters, respectively.
(a)(b)(c)(d)With everything considered, the results prove the viability of the cluster growth model for possible application on predicting the prime characteristics (Φ and I) of the occupied cluster in lattices (or in simple connected graphs) for the low-pressure range 0≤p<pc.


## 5. Conclusions

We considered a cluster growth model to predict the expected size of the invasion percolation clusters via an iterative calculation. A thorough description of the iteration was given. The method exhibits the phase transition of percolation models as the iteration diverges forp≥pc. The cluster growth model can only predict the size of the fully evolved cluster for 0≤p<pc as it cannot treat finite lattices yet. We provided a detailed interface update formula for the site square lattice and extensively tested the cluster growth model with this lattice type.We showed an application of the cluster growth model for the porosimetry percolation simulation. A significant finding of this study is the proof that the mapping between the input and the output of the simulation only depends on the topology of the network porosimetry percolation. This allows us to study the problem with uniform resistance distribution without the loss of generality.Calculations with the cluster growth model and porosimetry percolation simulations were carried out for different scenarios to compare the calculated expected cluster size with the simulated mean cluster size. We performed calculations and simulations assuming connected starting sets. Though the cluster growth model does not take into account the finite size of the lattice, the results still show that the evolution of the simulated and calculated cluster sizes nicely agree. The cluster growth model mostly provided an accurate estimation of the cluster size in the low-pressure range0≤p<pc. For instance, the mean relative errors of the iterative cluster growth calculation against the simulation results were 2.6−7.9% for L=1000,2000 lattices for the 1-sided and 4-sided invasion cases. The RMS of the relative errors were 3.4−11.4% for the same cases. A significant finding of this study is the proof that the mapping between the input and the output of the simulation only depends on the topology of the network porosimetry percolation is simulated on.In future work, we intend to apply this model to networks having different topology, e.g., three-dimensional lattices. To do this, an extensive number of numerical simulations will be required to determine an accurate relation between the cluster size and the interface size. We also plan to improve the cluster growth model to manage finite lattices. Then, we will calculateΦ and I for pressures p≥pc when a finite lattice is considered.The far aim of this work is to test and utilize the proposed cluster growth calculation on simple, connected graphs representing real rock samples with real pore network statistics.



---

*Source: 9919333-2021-07-09.xml*