NMR flow mapping and computational fluid dynamics in Ising-correlated percolation model objects

Water flow through quasi two-dimensional percolation model objects was studied with the aid of NMR velocity and acceleration mapping techniques. The model objects were fabricated based on computer generated templates of Ising-correlated percolation clusters of different growth/nucleation ratios for the occupation of the base lattice sites. The same pore networks were used for computational fluid dynamics simulations of hydrodynamic transport properties including hydrodynamic dispersion of tracer particles. The percolation threshold turned out to be lower than that in the uncorrelated case and adopts a minimum if cluster growth is about 100 times more likely than nucleation of the new clusters. The experimental and simulated flow velocity and acceleration maps coincide in great detail. The data have been analysed in terms of histograms and spatial autocorrelation functions. Furthermore, the travelling time of a tracer particle across percolation clusters was evaluated as a function of the Péclet number.


Introduction
The objective of this study is to model the pore network of porous media by 'Ising-correlated percolation clusters' in contrast to purely 'random percolation clusters'. This is an attempt to account for the growth history on which natural porous media are often based. Hydrodynamic transport of fluids through pore networks is of general interest. Apart from the model aspect, we therefore focus on experimental and computational techniques for the elucidation of characteristic transport features.
Modelling of interconnected pore spaces of porous media with respect to fluid transport properties raises many intriguing problems. Discussions can be found in [1]- [6]. Studies published in the literature focus mainly on random percolation clusters [5] and are often based on computer simulations [7]- [13]. Relatively little experimental work has been done so far in this respect. Before this background, the objective of the present study was (i) to examine Isingcorrelated percolation clusters generated with the aid of a Monte Carlo method, (ii) to elucidate fluid transport properties with the aid of a computational fluid dynamics (CFD) program and (iii) to compare the numerical results with experimental magnetic resonance imaging (MRI) techniques [14] applied to percolation model objects [15]- [21].
The problems encountered in comparisons of numerical simulations and experimental data in the context of percolation problems pertain (i) to finding appropriate evaluation protocols and characteristic parameters of maps of transport quantities and (ii) in having correlation lengths large enough so that the so-called scaling window applies where fractal power laws can be expected [7]. The latter point means that the porosity should be very close to the percolation threshold which in turn makes experiments difficult to perform.

3
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Nuclear magnetic resonance (NMR) offers a number of imaging variants permitting one to map transport quantities of interest quantitatively. Spin density maps can be used to display interdiffusion of isotopically labelled fluids selectively filled into adjacent compartments of porous model objects [19] (apart from rendering images of fluid-filled pore spaces in general). Velocity (vector) mapping and acceleration (vector) mapping visualize hydrodynamic flow arising from pressure gradients [15]- [18], [22]. Thermal convection flow can be examined likewise with the same sort of experimental protocol [21]. Temperature distributions due to thermal convection and conduction in porous materials can indirectly be probed by calibrating temperature-sensitive quantities such as the spin-lattice relaxation time accordingly [21]. Hydrodynamic dispersion coefficients characterizing the interplay between (incoherent) selfdiffusion and (coherent) flow can be studied using pulsed field gradient spin echo techniques [23]. A further method of interest is electric current density mapping of ionic currents in electrolyte solutions filled into a pore space in the presence of an electrostatic potential gradient [20]. Recent reviews on the perspectives of the magnetic resonance methodology for studies of porous media can be found in [24]- [26]. The next problem refers to the evaluation of experimental or simulated maps of transport quantities. Evaluated relations permitting the quantitative determination of characteristic parameters for the percolation properties are volume-averaged measurands (pore space density, velocity, acceleration, current density etc.) as a function of the porosity, histograms and spatial correlation functions. In particular, if the data reveal scaling law properties, a number of well-defined parameters can be evaluated from such secondary plots. Examples are the fractal dimension d f , the correlation length ξ and the percolation probability P ∞ characterizing the dependence of volume-averaged quantities on the porosity [5]. Correlation lengths can be evaluated from spatial correlation functions. Histograms permit one to identify different processes occurring on different scales of the measurand. If a description based on scaling laws applies, the determination of power law exponents can be of particular interest for theoretical considerations [4,5].
The strategy we pursue in the present study of fluid flow through correlated percolation clusters resembles the modus operandi of our previous investigations of random percolation clusters which were of 'random site', 'swiss cheese' and 'inverse swiss cheese' type [15]- [23]. The fabrication of well-defined model objects based on computer-generated templates provides not only the samples for informative NMR imaging experiments but also the complete datasets describing the pore spaces of these samples in all details. On the basis of these datasets, reliable CFD simulations using finite element method (FEM) or finite volume method (FVM) software packages can be carried out for comparison with experimental findings.
The advantages of this sort of system are obvious: in principle, the pore space is of sufficient complexity to represent natural or technical pore networks of interest. Since the model objects are based on computer-generated clusters, the pore spaces are well defined so that point by point datasets describing the pore space are available. Since these datasets are known, they can directly be fed into finite-element or finite-volume CFD programs in order to simulate transport properties. The percolation model objects are thus taken as a transport paradigm for any pore network of major complexity.

The Ising-correlated site-percolation model
The correlated site-percolation model considered in the present study is based on the Ising model in two dimensions [27]. A quadratic base lattice is assumed. Neighbouring sites form a pore if the quadratic pixels attributed to them share a common edge. The lattice sites are occupied with a random-number algorithm until the desired porosity is reached. The occupation probability of an unoccupied site depends on whether any of the neighbour sites is occupied or not. That is, nucleation (i.e. there are no occupied neighbour sites) is distinguished from growth (i.e. there is at least one occupied neighbour site) occupation probability. If both occupation probabilities are equal, the result will be a random site-percolation cluster. If the growth occupation probability is larger than the nucleation occupation probability, an 'Ising-correlated percolation cluster'results. The cluster-generating algorithm is specified in figure 1. The program is initiated by defining a square matrix of L × L lattice sites, the porosity to be finally reached, p final , the nucleation probability w nucl and the growth probability w growth . The pixel matrix is then occupied in an iterative way until the desired porosity is reached. In each iteration step, a pixel is selected randomly. If it is still unoccupied, the program continues to check for any occupied neighbour pixels. Depending on whether there is an occupied neighbour or not, the pixel is occupied with the growth or nucleation probability respectively. The ratio between the two occupation probabilities will be termed the 'growth-nucleation quotient', 'Black' means 'unoccupied', that is the matrix. White pixels represent the pore space. The base lattice size was 100 × 100 sites.
A growth-nucleation quotient f c = 1 means complete randomness of the pixel occupation. In the opposite limit, f c → ∞, the correlation of occupied pixels is complete so that just one cluster comprising all occupied pixels exists.
In figure 2, a typical example of an Ising-correlated percolation cluster is juxtaposed to a random site-percolation network generated for the same porosity with the aid of the algorithm given in figure 1. The concentration of the occupied pixels in compact clusters is characteristic for large growth-nucleation quotients.
The examination of Ising-correlated percolation clusters discussed in the following refers to four different square base lattices (L × L = 50 × 50; 100 × 100; 200 × 200; 500 × 500) in order to identify size-dependent effects. For each of these base lattices 100 independent clusters were generated. The parameters evaluated from these clusters were averaged accordingly. Different final porosities were investigated and the growth-nucleation quotient was varied in 10 steps in the range 3.2 f c 100 000. The total number of clusters considered in this study was 12 000 in various combinations of the cluster-defining parameters. Three of the clusters were examined in more details with the aid of CFD simulations and, after the fabrication of the corresponding model objects [15]- [21], with MRI experiments [14]. It is needless to say that these samples are based on 'infinite' clusters connecting two opposite sides of a square matrix. Figure 3 shows an example.

Percolation threshold
The percolation threshold [5] was determined as a function of the growth-nucleation quotient on the basis of 100 independent networks generated in each case for a given parameter triple L, p final , f c . The number of percolation networks containing an 'infinite' cluster, n 'infinite' , the so-called 'percolating' or 'sample spanning' cluster, was counted and related to the number of networks without percolating cluster, n finite . The ratio n finite /n 'infinite' was then plotted versus the porosity p final . An example is shown in figure 4. The percolation threshold is defined as the position where the ratio n finite /n 'infinite' takes the value 0.5. Figure 5 shows the threshold values as a function of the growth-nucleation quotient. The values tend to be lower than in the uncorrelated case (p uncorr c = 0. 592746 for random site percolation on a square base lattice; compare also [5,28,29]). Interestingly, a minimum appears in the range 30 < f c < 320. The minimum value is between 0.51 and 0.52 independent of the base matrix size as far as the investigation of this study is concerned. The literature value for the random site-percolation threshold in two dimensions, i.e. in the limit when f c approaches the value 1, is reproduced by the analysis used in the present study with an accuracy reasonable in the light of the small base lattice sizes to which our study was restricted. Note that finite percolation networks theoretically have thresholds somewhat larger than infinite ones [5].

Characteristic structure parameters of the percolating clusters
The correlation length has been evaluated from the volume-averaged porosity as a function of the probe volume radius [18]. For the determination of the volume-averaged porosity we have employed the so-called sandbox method [4,15]. The volume-averaged porosity is defined by where r | r k − r j |. The 'pore space density' ρ can take two values according to Equation (2) means that disk-shaped probe areas ('probe volumes') of radius r are centred in positions r k within the pore space. That is, r k must be the position of an 'occupied'pixel ('voxels'). N p (r) is the number of probe areas of a given radius r that can be placed with centres at pixels in the pore space. The probe areas must not exceed the sample under consideration, of course. N V is the number of pore space pixels within a given probe area. So, we first average the 'pore space density' ρ over the 'volume' of a given probe area of radius r (sum index j) and then take the arithmetic mean of the results for all possible probe areas of the same radius (sum index k). The volume-averaged porosity as a function of the probe area radius decays from the value 1 at a radius corresponding to a single lattice constant to a terminal value P ∞ , the 'percolation probability' at radii much larger than the correlation length, r ξ. For 'fractals', the initial decay can be described by a power law according to where d f is the fractal dimension, P ∞ is the percolation probability defined as the probability that a site belongs to the percolating cluster traversing the whole sample, ξ is the correlation length, b is the base lattice constant (≈0.4 mm for the present objects) and d E is the Euclidean dimension (= 2 in the present case). Double logarithmic plots of ρ V (r) can thus be evaluated with respect to the three characteristic parameters d f , P ∞ and ξ. Figure 6 shows the correlation length as a function of the growth-nucleation quotient for a relatively large system based on an 500 × 500 site lattice. Within the range of the examined growth-nucleation quotient, the data for a L × L site system can be described by the power law The parameters β ≈ 0.32 ± 0.04 and γ ≈ 7.75 ± 0.06 are constants. The correlation length 'increases' with the growth-nucleation quotient. There is no or only a negligible dependence of the correlation length on the porosity in the investigated range. The reason is that the number of clusters realized in the generation process depends on the nucleation probability but not on the porosity. If the nucleation probability is low, clusters once nucleated will preferably grow at the expense of the nucleation of new clusters, so that the correlation length becomes larger with lower nucleation probabilities.
On the other hand, the power law given in equation (5) is only valid as long as the correlation length is much shorter than the system size. Otherwise a 'cut-off' of the correlation length versus growth-nucleation quotient curve is observed. For clusters on a base lattice of size L × L = 50 × 50, for instance, we found a maximum correlation length of about 25b reached at a growth-nucleation quotient of about 3000. Beyond this growth-nucleation quotient value, the condition b ξ L is violated and the evaluated correlation length does not increase any further.

Fabrication of model objects
The templates for the model objects used for the MRI experiments were generated with the algorithm given in figure 1 using the IDL5.3 software package (Interactive Data Language based on a C++ platform). Since only the percolating (sample spanning) cluster is relevant for transport, all isolated clusters were removed (compare figure 3(a)). The porosity is nevertheless specified by the value p final anticipated in the computer generation process (see figure 1). A comparison of the porosities p final of the whole percolation networks with the porosities effective for the sample-spanning cluster alone, p cluster is given in table 1.
The datasets resulting from the cluster-generation program, figure 1, were fed into a digital circuit board plotter (LPKF, 'Leiterplatten-Konturfräsen') used for milling the pore space into 3 mm thick polystyrene plates. The milling depth was 2 mm and the location resolution of the milling tool was 6 µm (see figure (3b)). Porous model objects have been fabricated for growthnucleation quotients 10, 100 and 500. The base lattice size is L × L = 200 × 200 sites in all cases. The resulting porosities are p final = 0.59, 0.57 and 0.60 respectively.  The final, quasi-two-dimensional model objects are composed of stacks of nine identical plates structured in this way in order to improve the signal-to-noise ratio of the NMR experiments. The stacks have been incorporated in plexiglass containers with liquid inputs and outputs (see figure 7). The pore space was filled with demineralized water doped with CuSO 4 in order to reduce the spin-lattice relaxation time T 1 to about 500 ms. A typical spin density map recorded with one of the samples is shown in figure 3(c). Note that the matrix material and the container do not provide liquid-state NMR signals and can therefore easily be discriminated from the liquid.

NMR imaging system
Spin density, velocity and acceleration mapping experiments were performed with a home made NMR imaging system using a Bruker 4.7 T magnet with a 40 cm horizontal room temperature bore. All NMR signals refer to proton resonance at 200 MHz. The probe head was a Bruker birdcage resonator with an inner diameter of 195 mm. The maximum field gradient

velocity (trains a plus b) and acceleration (trains a plus c) mapping experiments with quasi-two-dimensional objects in the x-y plane. A standard
Hahn spin echo imaging sequence is used in combination with one (b) or two (c) pairs of gradient pulses (see [14,21]). The arrows indicate the mutual polarity of the gradients which are incremented in nine steps in a corresponding series of transients. The ramps of the gradient pulses are not shown for simplicity. The sliceselective π pulse ensures that no signals from the fringe areas above and below the objects are recorded. Instead of the positive compensation gradient pulse for frequency encoding along the z-direction, a negative gradient lobe (dotted line) just before the read-out gradient may be preferable in order to avoid motion artifacts in the imaging part of the pulse sequences. To acquire the velocity or acceleration vector field, all components must be probed one by one. For quasi-two-dimensional objects of the present investigation, the dataset acquired in this way consists of two four-dimensional matrices. Fourier transforms in all four dimensions lead to conjugated datasets providing the two in-plane velocity or acceleration vector components. In the velocity or acceleration maps, the magnitudes v = v 2 x + v 2 y or a = a 2 x + a 2 y are plotted where the flow directions are perpendicular to the z-direction.
was 50 mT m −1 in all three space directions. The experimental spin density maps reproduce the computer-generated template structure reliably (see figure 3) although the quasi-two-dimensional model objects are subject to friction at the bottom and top surfaces whereas the simulations refer to two dimensions in the proper sense. This coincidence demonstrates the quality of both the fabrication and NMR imaging processes.
Spin density maps of the water-filled pore spaces were recorded with the aid of the radio frequency (RF) and field-gradient pulse sequence shown in figure 8(a) [14]. The maps were rendered by processing the raw data by two-dimensional Fourier transformation. Due to field-of-view limitations, the largest objects that can be mapped with our system are restricted to a base lattice size of 200 × 200 or an object size of 10 × 10 cm 2 . For the flow experiments, a stationary pressure gradient was exerted on the fluid-filled pore spaces of the model objects by connecting a reservoir 1.5 m above the sample level. The reservoir was permanently refilled with the aid of a peristaltic pump. The total flow rate was in the range 0.1-0.3 ml s −1 .
Phase encoding of flow velocity or acceleration components was achieved by applying additionally the field gradient pulse sequences shown in parts (b) and (c) of figure 8, respectively. Detailed descriptions of velocity and acceleration NMR measuring techniques can be found in [14,21,30]. The digital space resolution was better than 300 µm in all three space directions. All experiments were carried out at room temperature. Typical velocity and acceleration maps recorded with the object shown in figures 3(b) or 7 are rendered in figures 9 and 10(c).

CFD software
All flow experiments were accompanied by CFD simulations based on the same pore space structures as in the experiments. The simulations were performed on a PC using the software

Nuclear magnetic resonance imaging experiments
The experiments were performed under conditions of locally stationary velocities and accelerations. Velocity and acceleration vectors of the spin bearing particles were recorded at certain positions (voxels) rather than for a certain tracer particle experiencing varying velocities and accelerations when travelling through different voxels. Stationary flow patterns and pressure gradients stipulate locally stationary velocities and accelerations. Figures 9 and 10 show examples of experimental maps of the velocity or acceleration magnitudes. The flow patterns are reasonably well reproduced by the CFD simulations displayed in figure 10. Note the discontinuous distributions of the acceleration along the flow pathways (figure 9, right column), whereas the velocity is represented by continuous patterns (figure 9, left column).

Simulations
Maps of the velocity magnitude, v = v 2 x + v 2 y , in the two-dimensional percolation clusters were obtained as described previously [18] by solving the conservation equations for mass and momentum (see [31]) on a discrete domain. Examples are shown in figure 10(a).
Hydrodynamic dispersion of a neutral tracer particle diffusing in the streaming fluid was simulated by taking the equation of energy into account. Typical trajectories are shown in figure 11. The transport of the tracer depends on the macroscopic Péclet number [2] defined by where v is the fluid velocity, D is the molecular diffusivity of the tracer and l is a characteristic length of the flow field. In site-percolation models, l can be identified with the correlation length ξ.

Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
That is Here we anticipate that the flow velocity is constant over the correlation length for simplicity. In the absence of flow, the tracer executes a random walk of step size s = b/10, where b represents the minimum dimension of the pore size which is 400 µm in our case. The tracer trajectory also depends on where it is launched into the flowing fluid. An example based on different Péclet numbers is shown in figure 11. For the lowest Péclet number, that is Pe = 0.09, diffusion strongly determines the tracer trajectory. On the other hand, if the Péclet number takes a high value such as Pe = 40.35, flow convection dominates the dispersion almost completely.
Maps of the acceleration magnitude, a = a 2 x + a 2 y , cannot be simulated directly with the FLUENT 5.5.14 software package. However, velocity maps obtained in the usual way can be processed correspondingly using the IDL 5.3 software mentioned above. Anticipating an Eulerian approach, the two-dimensional velocity vector field v = (v x (x, y, t), v y (x, y, t)) can be differentiated according to Note that the acceleration components given in the equations (8) and (9) depend on the (voxel) position but not on time since all velocity components are assumed to be locally stationary and merely change from position to position. Examples are shown in figure 10(b).

Evaluations of the experimental and simulated maps
Maps of the (water) spin density, velocity magnitude and acceleration magnitude were rendered for the three examples of Ising-correlated percolation clusters as described above both with NMR experiments and computer simulations. In each case, the experimental and simulated patterns reproduce each other reliably. Figures 9 and 10 show typical velocity and acceleration maps obtained from NMR experiments and CFD simulations respectively. The velocity maps indicate pathways of a more continuous character, whereas the acceleration strongly fluctuates along the flow channels. That is, the acceleration is only in certain voxels large enough to be detected above the noise level of the experimental set-up. These sections correspond to bottlenecks and strongly curved pathways where flow is strongly accelerated or decelerated.

Volume-averaged porosity
The volume-averaged porosity defined by equation (2) was evaluated from water spin density maps. An example is shown in figure 3(c). Evaluations of the fractal dimension according to equation (4) lead to values of d f ≈ 1.75 without significant variation for growth-nucleation quotients in the range 10 < f c < 1000. This value differs from the fractal dimension of uncorrelated clusters of the same porosity d f ≈ 1.86 evaluated from simulations or d f ≈ 1.83 evaluated from experiments by less than 6%, that is, within the error limits. Similar evaluations for volume-averaged velocity and volume-averaged acceleration defined analogous to equation (2) were suggested in [18]. A power-law decay to a long-distance plateau similar to equation (4) can however be stated only for low correlations corresponding to f c < 10. That is, this sort of behaviour appears to be valid only for clusters approaching the entirely random limit.

Transport backbone
Transport backbones can be evaluated from the maps by blackening all voxels with velocities below the noise level. The latter was determined in stagnant zones where there is no flow at all. The procedure is described in [15]. Typical maps of backbones are shown in figures 9 and 10. Taking the corresponding pixel subset of the spin-density maps, one can evaluate the fractal dimension of the backbone. The value, d bb f ≈ 1.25, is much smaller than that of the whole percolation cluster [18] as it must be since the pore space contributing to transport is only a small fraction of the total percolation cluster. evaluation error. The characteristic features reported previously for convection in uncorrelated percolation clusters [21] show closely up again and appear to be scarcely influenced by the growth-nucleation quotient. Similar to the cases described in [21], an exponential variation of the velocity distribution can be stated. However, unlike the thermal convection case in uncorrelated clusters, there is now some weak tendency towards bi-or multi-exponential decays (see figure 12(a), for instance) which is assumed to reflect the wide range in the size of the flow obstacles caused by the Ising-correlation.

Mean tracer travelling time
as expected for the limit where displacements by hydrodynamic flow dominate and where an effectively constant displacement rate along the net transport direction occurs on a length scale beyond the correlation length of the system. Below Pe ≈ 1, the steep decay of the mean travelling time with the Péclet number indicates the decreasing dominance of Brownian motions and the increasing influence of coherent flow.

Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
The spatial velocity autocorrelation function can be described in a wide section of this regime by a power law which is the more pronounced the larger the correlation length, that is, for larger growthnucleation quotients. The exponent α increases with the growth-nucleation quotient f c . A reliable evaluation of the exponent is however difficult for small f c values because region II then becomes rather narrow. The decay in region II gets steeper with increasing f c . This is considered to be a consequence of the size of the obstacles around which the liquid needs to flow. The obstacle size increases with f c (for an illustration see figure 2). The flow tortuosity caused by obstacles grows with the obstacle size at constant porosity. The velocity direction consequently varies in a wider range for larger f c whereas flow streamlines in entirely random clusters of high porosity tend to be straight with little directional fluctuations.
The 'plateau region' III arises for ξ r l cluster . Actually, the plateau character becomes more pronounced the larger the width of this region, that is for low f c . Otherwise finite-size cutoff effects and the influence of inflow and outflow at the edges of the cluster become perceptible and affect the plateau character of this region.

Discussion and conclusions
The above results refer to a combined experimental and computational study of the structure and transport properties of Ising-correlated quasi-two-dimensional percolation model objects. It is shown that characteristic parameters such as the percolation threshold, correlation length and fractal dimension deviate from the values found for uncorrelated, that is random site-percolation model objects reported in previous papers. Dependences on the growth-nucleation quotient and on the porosity have been examined.
Taking into account the finite size of the systems under examination here, the percolation threshold value is found to be lower than in random site-percolation networks and interestingly shows a minimum in a range 30 < f c < 320 of the growth-nucleation quotient. This is in accordance with findings reported in [9,28] for other correlated percolation models. Correlations significantly diminish the percolation threshold as a result of the clustering of, on the one hand, low and, on the other hand, high permeability regions. This is consistent with experimental observations that indicate a high permeability of many natural porous media even at very low porosities [2].
The explanation of the minimum of the percolation threshold is that, with increasing growth-nucleation quotient, the cluster size increases which favours transport across the network. However, there is an opposite tendency by the fact that the number of clusters is reduced at the same time. At growth-nucleation quotients far below 100, the first tendency dominates whereas above that value the second trend starts to govern the dependence on the growth-nucleation quotient.
The correlation length obeys a power law, equation (5), in a wide range leading to strongly increasing values for larger growth-nucleation quotients. This tendency is also visible in region II of the spatial velocity correlation function: This 'scaling region' is getting wider with increasing growth-nucleation quotient and is terminated by a crossover to a weakly indicated 'plateau region' at the correlation length (see figure 14).
The hydrodynamic dispersion of a neutral tracer in the flow field was studied as a function of the Péclet number using again computer simulation tools. Two regimes of the mean travelling time could be identified corresponding to an isotropic dispersion range, Pe < 1, where incoherent Brownian motions dominate, and a mechanical dispersion interval, 1 < Pe < 10 4 , where coherent flow controls dispersion almost completely [2,10].
While random site-percolation clusters show some universality with respect to the choice of the elementary pore shape ('random site', 'swiss cheese' and 'inverse swiss cheese') [18], the Ising-correlated structures considered here appear to form another category with different parameter features. Ising-correlated percolation clusters are expected to mimic the structure of natural porous materials that were formed by growth processes of a similar nature. The main pore space feature that matters here is the larger pore connectivity resulting from the correlation assumed in the growth algorithm generating the models. This is in contrast to random sitepercolation networks which are known to be too 'random' for the description of typical porous media.
As a general remark, it may be noteworthy in this context that the Ising-correlated quasi-twodimensional percolation networks considered here can be associated to the growth dynamics of urban clusters. Corresponding analyses have been published in [32,33]. This sort of comparison may be of interest with respect to the dependence of the cluster parameters on the growthnucleation quotient.