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 threshold . We describe an application for the porosimetry percolation model. The calculations of the cluster growth model compare well with simulation results.
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 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 , 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 and the type of the lattice [4–6]. It was found that there is a strict critical percolation threshold characteristic of the lattice type. For an infinite lattice, the probability of the existence of the spanning cluster is zero for , while it is one for . At , there is a singularity in percolation and many properties of the largest cluster follow a power law of . To this date, was derived or computed for many different lattice types. For a few special lattices, the exact value of 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 is solved exactly [9]. For other lattice types, an extensive number of simulations were carried out to determine numerically [6, 10].
The existence of the critical threshold 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 , 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 and has “immunity” otherwise. Similar to the concept of herd immunity, we can contain the finite cluster sizes, if [14].
We focus on the site variant of the problem, since it is more general than the bond variant (see Section 4.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 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 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 pressure (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 will be occupied. A site is accessible only if it is part of a path of interconnected sites in which every site has and this path is connected to the starting set.
In porosimetry percolation, we can track the evolution of the cluster as is increased, and compute, e.g., the cluster size as the function of . 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 . 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 . Can we also use the cluster growth model to calculate this expected cluster size for any ? 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 .
This article is structured as follows. In Section 2, 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 set of occupied sites in , also called a lattice animal or polyomino. The size of the polyomino is the cluster size .
Another important characteristic is the interface of the cluster. The interface is the boundary of , consisting of the perimeter sites, i.e., the set of sites that are not in but adjacent to some sites in . The perimeter is the number of interfacial sites (or perimeter sites), i.e., .
Let denote the number of lattice animals with size and perimeter . The so-called perimeter polynomial is the generating function:
For example, the first few perimeter polynomials for fixed polyominoes are as follows [24]:
The perimeter polynomials up to are published on Mertens’ webpage [25]. The number of lattice animals of size on the site square lattice is
The relation between the cluster size and the perimeter 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 is
As the enumeration of is still an open problem [26, 29], 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]:
For the mean interface size , there is an empirical formula for large polyominoes proposed by Conway and Guttman [33]:
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 ) for a prescribed pressure such that any site with cannot be occupied.
The probability that a cluster of size develops from a single starting site is given by the following [34, 35]:
This equation is the link between percolation clusters and lattice animals. Based on (8), one could derive (and ) as the function of 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 . Though invasion percolation is different from ordinary percolation, it also exhibits the same criticality feature. This means that above the critical pressure , that is equivalent to the percolation threshold, of the invasion becomes unbounded. For , the occupied cluster consists of a finite number of sites, but for , it becomes infinitely large. In particular, the percolation threshold for the site square lattice is [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 ’s sufficiently close to ; that is,
The “low-density” (small ) series expansion of the expected cluster size as a function of is in the following form (see, for example, [24]):where are constant coefficients. We determined these constants up to based on the perimeter polynomials.
In Figure 1, we depict the expected cluster size as a function of for the site square lattice. For the simulations, a finite lattice was created, and the starting set was an innermost site.
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 ’s but diverges from the simulated curve around . This is why we also calculated the Padé approximant [38] of order of (10) because it is accurate on a broad pressure range.
There also exists a formula similar to (9) for the mean interface size of large clusters and ’s sufficiently close to that was derived by Stauffer [34] and Wolff and Stauffer [39]. They suggested the following relation:where and , andis the scaling function [39]. We can also compute and from cluster size and interface size statistics of percolation clusters. We carried out 235000 numerical simulations on site square lattices and extracted the mean interface size 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 . Figure 2 shows the mean interface sizes corresponding to different cluster size intervals that are represented by their mean cluster size.
(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 and , there is also a softening nonlinear component that is noticeable at smaller ’s, see Figure 2(b).
Figure 3 demonstrates the difference between the perimeter distribution of all possible polyominoes and simulated percolation clusters of 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 pressure , the starting cluster size and interface size of the prescribed starting set, and the function . We want to calculate the evolution of the expected cluster size and the corresponding mean interface size . 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 Section 4.1) and a connected starting set. The expected number of sites that get occupied from the initial interface is at pressure ; this leaves interfacial sites unoccupied. After the first round of iteration, is obtained immediately as
We have 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 . In the next step, the cluster can only occupy the new interfacial sites; thus, the expected number of sites annexed by the cluster is . Figure 4 explains how the cluster and its interface evolve in one iteration and how to interpret .
(a)
(b)
Let us assume that and are already obtained as in Figure 4(b). The following sites constitute :(i)All the “extinct” interfacial sites that already contributed to (dots in Figure 4(b))(ii)“New” interfacial sites that became part of the interface in the th iteration (circles in Figure 4(b))
The extinct interface sites were already tested in the th step. Since they were not occupied, their resistance must be larger than . Therefore, in the th step, only the new interface sites will be considered for possible occupation. Based on Figure 4(b), can be derived and it is
Thus, the expected number of sites annexed in the th iteration step is . Now, we have and from this train of thought, a simple recursive formula is derived for the expected cluster size and interface size:where is a function of the cluster size. Based on the preliminaries (see Section 2), we expect that consists of a linear and a nonlinear part, i.e.,
We also expect to be a strictly monotonously increasing function. In Section 3.2, we dissect further for the site square lattice. Finally, the complete formula of the iteration describing the cluster growth model is as follows:
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 that holds); i.e.
Equation (18) can be transformed into a set of two first-order nonlinear difference equations as follows:
Let us investigate the stability of (19) via fixed-point analysis. Linearization of (19) around the fixed point of the iteration yields
The eigenvalues of are obtained by solving the following characteristic equation:
Based on Figure 2, we expect a softening nonlinearity; i.e., we have
We get the eigenvalues , for . Since , the fixed point is a neutral point for (Lyapunov stable). This latter can only hold for , so this is a condition for marginal stability. For , the second eigenvalue regardless of the value of and the system becomes unstable. This shows that the radius of convergence of the system (20) is for large ’s. The cluster growth model is consistent with the existence of the critical threshold , ifholds near . This result is also consistent with equation (11).
3.2. Interface Update for the Site Square Lattice
We attempt to construct a function 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 size of the cluster of expected size can be calculated with Stauffer’s formula (11). Substituting equation (12) into and then into (11) leads to
This equation also shows softening nonlinearity as the exponents . We recall that this formula is only accurate for large ’s near . To make use of (24), we have to determine the lower threshold 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 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 ; 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 against .
Figure 5 shows that is not decreasing for some ! Via iterative search, we found that has a maximum at (marked by a vertical line in the figure) and the corresponding threshold cluster size is . This shows that (24) is only applicable for .
Moreover, note that, in the limit of , equation (24) reduces to
We see that the applicability of (24) is heavily restricted: it is only valid for a single-site starting set, and . Moreover, we cannot use it, until the expected cluster size exceeds the threshold size . However, we need to update also for smaller ’s and for different starting sets in the iteration.
We assume that below the threshold , the cluster develops through separate branches. The threshold value corresponds to the particular starting set , with . We adjust for any larger values of , , and . A larger must increase the threshold size by the same amount, while a larger and/or increases the number of branches emanating from the starting set by the ratio . Equation (26) assigns accordingly:
Let us assume that the starting set is small, i.e., , but it is not restricted to the single-site case with , . To calculate and with high accuracy, a piecewise is proposed. The linear part is
The nonlinear part is
Explanations of the cases are as follows:(i)For , the expected cluster size is between that of a size of a single site () and a domino (). In this case, the interface size is trivial; see equation (6): .(ii)In the range , we use a fit obtained from the mean interface size vs. cluster size statistics from percolation simulations; see Figure 2. The constant parameter 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 , we can finally use Stauffer’s formula (24). The constant parameter 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 . For such large starting sets, we will simply use
That is consistent with the limit expression (25) of Stauffer’s formula for large cluster sizes. The relaxation of is required at some iteration steps to avoid two possible contradictions. For starting sets with , it is possible to get for low ’s. This is highly atypical during the low-pressure build-up of percolation clusters. Moreover, we experienced that Stauffer’s formula may yield when , whereas ratio must decrease as increases; see Figure 5. Therefore, we also incorporate a simple relaxation:
Now, we show an example of how and evolve in the cluster growth model. We study the plane depicted in Figure 6. All the possible 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 for pressure and starting . By design, the iteration attracts any pair of initial values onto in one iteration step. Thus, the iteration jumps to located on and then follows towards the fixed point . Upon reaching , the function 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.
4. Application for Porosimetry Percolation
A porosimetry percolation simulation can compute the evolution of the cluster for any , can be computed for any , 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 , since the iteration diverges for . 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 , the saturation curve shows as a function of . 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.,where is the side length of a finite site square lattice. The Washburn equation [40] specifies the relation between pressure and pore size of the smallest pore, which can be filled up with mercury at that pressure. That is, a larger is required to fill a pore characterized by a smaller . The pore size is a physical property that is interpreted in percolation models as , the resistance of the pore. We clarify that resistance is a pressure-type quantity; it is equal to pressure 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 . As an initial guess, the PSD derived from 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 . The main idea was to treat the simulation as an input-output mapping: from input , it yielded output . The task is to find input for which the output is the measured . Then, the resistance distribution and the PSD can be derived from instead of . The method is illustrated in Figure 7.
The mapping relates and as
Since the saturation curves are strictly monotonously increasing functions on , is a continuous bijection, whose inverse is also continuous; i.e., is a homeomorphism. Thus, the inverse of exists and it can be used to obtain as
It was argued in [20] that does not depend on the input , 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 and the starting set is a single site. Since resistance is a pressure-type quantity, we can describe the input resistance distribution with a cumulative distribution function (CDF) .
Distribution shows the probability that the resistance of a site is . Based on equation (8), the probability that the occupied cluster will have a size of is given by
The expected cluster size is
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 sites is as follows:
We can calculate (36) for all to get the saturation curve .
Let us draw the resistances of the sites from a uniform distribution; i.e., we have . In this case, the identity function is transformed into by the mapping ; i.e.,
Do we get the same for any resistance distribution ? Based on (36) and (37), we can formulate
Substituting for into (38) yields
We got back equation (32) on the left-hand side with different notations. This proves that mapping is indeed independent of .
The independence of from 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 , e.g., using the proposed cluster growth model.
4.2. Tests with Small Starting Sets ()
We test the cluster growth model on the site square lattice with the interface update proposed in Section 3.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 for for three different starting sets: single site (), maximum perimeter tromino (), and maximum perimeter heptomino (). We also carried out porosimetry percolation simulations with these starting sets on a site square lattice of size (the number of sites was ). 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 vs. the iteration number is depicted for the single-site starting set; the parameter of the curves is .
(a)
(b)
Figure 8 shows that the evolution calculated with the cluster growth model is close to the simulated evolution for . Results are also shown for some cases to demonstrate that the iteration indeed diverges for those pressures.
In Figure 9, is shown for the pressure range 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 iterations, i.e., . The simulation results are averaged over 1000 simulations.
(a)
(b)
(c)
Figures 9 and 10 show that the cluster growth model provides a good approximation of for , though the calculation overshoots both the simulation result and Sykes’ expected cluster size close to . 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 for , , where were chosen. The mean and root mean squared (RMS) of these relative errors are listed in Table 1.
According to the relative errors, the cluster growth model considerably overshoots the simulation result near . Still, the relative errors are acceptable if we exclude the data points that are very close to . 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 , if the simulations were carried out on a much larger lattice.
In Figure 11, representative visualizations of the evolution on the plane are shown for . The evolution of the cluster growth model shows nice agreement with the simulated evolutions.
(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 range .
4.3. Results for Large Starting Sets ()
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 for inputs 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 , we have
In the 4-sided invasion case, the starting set completely surrounds the rest of the lattice; we have
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 (). In Figure 12, the expected saturation of the lattice is shown in the pressure range for the different starting sets. The calculations of the cluster growth model were performed with iterations, i.e., . The simulation results are averages of 1000, 500, and 100 simulations for , respectively.
(a)
(b)
(c)
(d)
(e)
(f)
The plots are convincing regarding the accuracy of the cluster growth model. Similar to Section 4.2, the relative error of the iterative cluster growth calculation against the simulation data was computed for , , where were chosen. The mean and RMS of these relative errors are listed in Table 2.
The 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 smallest case, the expected cluster size significantly overshoots the mean simulation result near . 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 Figure 13, representative visualizations of the evolution on the plane are shown for .
(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 ) of the occupied cluster in lattices (or in simple connected graphs) for the low-pressure range .
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 for . The cluster growth model can only predict the size of the fully evolved cluster for 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 range . For instance, the mean relative errors of the iterative cluster growth calculation against the simulation results were % for lattices for the 1-sided and 4-sided invasion cases. The RMS of the relative errors were % 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 for pressures 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.
Data Availability
All datasets generated for this study are made available on request through the corresponding author of this article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the NRDI Funds (TKP2020 IES, grant nos. BME-IE-WAT, TKP2020 NC, and BME-NCS) based on the charter of bolster issued by the NRDI Office under the auspices of the Ministry of Innovation and Technology.