© E. Shahrian & M. Masihi, published by IFP Energies nouvelles, 2021

Nomenclature

P_i : Initial pressure (psia)

_Pwf : Bottom hole flow pressure (psia)

Φ: Porosity fraction

r_w : Wellbore radius (ft)

μ : Fluid viscosity (cp)

K : Average permeability of the reservoir (mD)

q : Flow rate (STB/D)

B : Formation volume factor (bbl/STB)

c_t : Total compressibility (psia⁻¹)

m : Slope of the line in log chart

t : Time (hr)

h : Reservoir height (ft)

s : Skin

P(DD): Pressure at Datum Depth (psia)

DD: Datum Depth (ft)

P_b : Bubble-point pressure (psia)

Bo: Oil formation volume factor (bbl/STB)

γ_o : Specific gravity of oil

γ_w : Specific gravity of water

R_s : Gas solubility (scf/STB)

c_R : Rock compressibility (psia⁻¹)

T : Reservoir temperature (F⁰)

P′: Derivative pressure (psia)

p_c : Percolation threshold for an infinite system

: Percolation threshold of an infinite size system

β : Connectivity exponent

ξ : The correlation length of a system.

1 Introduction

Well testing is a popular technique for the evaluation of reservoir parameters such as reservoir permeability. Therefore, several well testing methods are proposed such as Drawdown Test, Buildup Test, Interference Test, Drill Stem Test (DST), Falloff Test, Multi-Rate Test, Multiple-well Tests, and Pulse Tests [1, 2]. However, among these the Drawdown test is a common and simple well testing method with significant results [3].

To do this, a shut-in well opens to produce the reservoir fluid with a constant flow rate. Herein, several assumptions are made such as single-phase flow, and constant fluid compressibility and viscosity, while the capillary pressure effect is assumed negligible. Often, the Drawdown test applies in the exploratory wells as well as newly drilled ones. It aims to evaluate the boundaries of the reservoir, the effects of the skin, the average permeability and pore volume of the reservoir. The advantage is that there is no loss of cash flow due to the shut-in as is the case with the buildup test. But, the downside could be the challenge of keeping the flow rate constant [4]. Significantly, the drawdown test is capable of detecting the reservoir heterogeneities including the existing boundaries around the wellbore or different layers with different permeabilities. This characteristic is among the most valuable features of this study.

Equation (1) shows the resulting pressure equation for the drawdown test in a wellbore located in a cylindrical reservoir with radial flow geometry and infinite-acting region [4]:

(1)

As is known, the slope of this equation on semi log scale gives the permeability of the reservoir as given in equation (2):

(2)

Additionally, the skin could be computed through equation (3):

(3)

1.1 Detecting reservoir layers and boundaries

During the test, if the reservoir size is small or the test period is long enough, then the radius of investigation reaches the boundaries of the reservoir. In such a case, the boundary effects would be observed in the well testing analysis and pressure signals could be helpful to estimate some parameters such as distance to the barriers or boundaries, and the shape of the reservoir. Quantitatively, the distance to the reservoir boundaries, the distance and strength of a close aquifer, the distance to a gas cap and the advancing waterfront, and finally the shape and size of sand bodies in the reservoir could be computed. On other side, this data could provide qualitative details on the nature of the reservoir and the flow behavior close to the wellbore.

As generally applied, two types of reservoir boundaries are: (a) impermeable and (b) constant pressure boundaries [5]. The impermeable boundary (also known as a closed boundary) is associated with volumetric reservoir while, a constant pressure boundary assumption produces an effect that is approximately close to the constant pressure limit [5]. The lapsed required time to observe the boundary effects mostly depends on the distance to the boundary, the permeable formation properties and the type of fluids.

The pressure response in a wellbore near a linear boundary could be analyzed using the superposition technique, which is also called the method of images. This analytical method uses production or injection wells to account for the effects of no flow or constant pressure boundaries. However, a more general approach to obtain the pressure response, especially for heterogenous reservoir model, is through numerical methods.

But the case is a little bit different for faults. The detection of faults in the reservoir is usually performed through geological and seismic studies. The fault itself may be sealing or leaky. The sealing fault completely obstructs the lateral flow of the fluids and may work as a part of the hydrocarbon accumulation. A single finite sealing fault acts as a no flow boundary type.

To detect the heterogeneity we need to differentiate between various heterogeneous reservoirs. Due to the sedimentation process, most oil and gas reservoirs were stratified to some degree. These reservoirs are usually divided into two groups: (1) layered reservoirs without crossflow (no vertical flow), where layers communicate only through the wellbore, and (2) layered reservoirs with crossflow between layers, where layers communicate at their contact planes throughout the reservoir [4]. Not to mention that a reservoir layer is considered a reservoir rock with homogeneous lithological properties on a macroscopic scale. Hence, a layered reservoir is divided into laterally continuous layers with differing permeabilities and porosities. For a non-communicating layered reservoir, the pressure behaves on well tests as if it is homogenous and isotropic. However, the behavior of the layered reservoir with crossflow is different from the behavior of a homogeneous and isotropic reservoir. The test objectives would be to determine permeability k and skin s for each layer and the average reservoir pressure p.

The conventional drawdown and buildup tests only are useless here as they only show the pressure behavior of the whole reservoir environment around the wellbore. Hence, the behavior of a multilayer reservoir may not be distinguishable from the behavior of a single layer formation, even though a multilayer reservoir may have a distinct behavior without wellbore storage effects [6].

But the pressure-derivative behaviors have shown of great applicability. As is illustrated in Figure 1, for the homogenous reservoirs, the derivative curve behaves as a horizontal straight line in the transient flow regime followed by a sharp decrease for the case of a constant pressure boundary, and a significant increase for the no flow boundary case. But, for a dual-porosity reservoir, a valley in the derivative curve is observed and for a triple-porosity model, two continuous valleys are recorded [7].

Fig. 1 Diagnostic plot of reservoirs. (a) Homogenous reservoir with constant boundary pressure, (b) homogenous reservoir with no-flow boundary, (c) dual porosity of the naturally fractured reservoir, (d) triple porosity of the naturally fractured reservoir [8].

2 Well testing in heterogeneous reservoir

The heterogeneity of permeability significantly impacts the fluid flow dynamics [9]. Neglect of this important parameter may lead to failure in an expensive EOR process or to leakage of a large volume of stored CO₂ [10]. Among different properties, permeability directly controls the flow behavior. Moreover, its variation is usually much greater than other property variations. These reasons explain why most of the heterogeneity measures concern the degree of permeability variation [11, 12].

Permeability distribution in a heterogeneous reservoir is a complex problem with coefficient of variation (i.e. the ratio of mean to standard deviation of permeability distribution) around or more than 1. The quality of its characterization depends on the number and types of data available including core samples and well test data. There are also other sources such as well logging, Geophysical data, and Geological Interpretations that could help to have a better estimation [13].

Heterogeneity classification criteria could be classified into two main categories: static and dynamic. The most famous static criterion is the Dykstra-Parsons coefficient of heterogeneity [11] that might be used for fluid flow studies including dispersion in porous media [14], recovery of oil reservoirs [15, 16], and carbon dioxide storage in aquifers [17]. The fluid flow information in dynamic criteria is also used to express the level of the reservoir heterogeneity. One such dynamic measure is dispersivity which depends on the fluid flow physics. Among the dynamic heterogeneity criteria are the Lorenz dynamic coefficient [18] and the vorticity [19].

One practical application of this phenomenon is on the well testing response. Gautier and Noetinger [20] proposed a method based on the Bayesian inversion theory, in conjunction with the fast evaluation technique of well testing, to estimate the geostatistical parameters such as the correlation length and the permeability variance from well test data. But, the effective permeability provided by well tests is often different from the “static” effective permeability. This is due to the difference in dynamic condition and the sample volume. To circumvent this problem, Noetinger and Gautier statistically derived the ensemble-averaged pressure in the medium from which the effective conductivity may be inferred [21].

As emphasized by Hewett (1986) the transport properties are mainly controlled by the details of spatial correlations in the permeability map. Al-Khidir et al. (2011) used fractal geometries to support the reservoir heterogeneity model [22]. Fractal geometry provides an appropriate mathematical model of the distribution of rock properties which uses a power law relation. The fractal based structure creates distributions that show both heterogeneity and ordering pattern [23].

Moreover, Sagar et al. [24] used analytical solutions for radially heterogeneous reservoirs to define an equivalent radial permeability for the hetrogeneoug region of investigation. This was compared with numerical experimentation on drawdown well-test simulations and concluded that the instantaneous well-test permeability can be explicitly related to the distribution of small-scale permeabilities within the annular regions of the reservoir.

Sureshjani et al. [25] studied the permeability variation from the analysis of pressure transient on non-radial flow regimes and proposed the use of pressure drop at the boundary of the region of investigation to account for the heterogeneity.

Reservoir heterogeneity can be seen on various scales from milli meters to kilometers. Upscaling approaches or scale up methods are aimed to cope with the impact of heterogeneity. Some averaging methods can measure small-scale heterogeneity (i.e., due to grain size), while large-scale heterogeneity (i.e., due to geological layers) is considered with appropriate reservoir layering. The most challenging parts in heterogeneity modeling are small to medium scale heterogeneity, such as low-permeability barriers (Lasseter et al., 1986). Although the knowledge of permeability distribution at smaller scales can be used to estimate an upscaled reservoir property (Hunt, 1998), the value is most likely related to the critical path (percolating) resistance [26].

The issue of heterogeneity and upscaling has also been studied by Matheron [27], Cardwell and Parsons [28], and Dagan [29] for 2D systems. He has reported that the effective permeability is some numbers between the harmonic and the arithmetic mean i.e. . As emphasized in Renard and De Marsily [30] the power-law averaging may give the effective permeability i.e. with power exponent p bring between −1 and 1. Noetinger found that a power relation with exponent works for statistically homogeneous and isotropic media [31]. The power value of p numbers depends on the spatial distribution of permeability. In 2D, it corresponds to the geometric average that was found to be exact in 2D for lognormal media by Matheron (1967), who derived it using an elegant duality argument specific to 2 dimensions [32].

In this work, the effect of these small-scale heterogeneities on the reservoir’s dynamic performance is studied. In particular, we concentrate on the pressure response from a production well located in a percolation based porous medium. But, we initially need to make such a heterogeneous porous medium and compare the pressure response results with the analytical/numerical solution of a similar problem on a homogenous porous medium. The next section covers this subject.

3 Building a percolation-based reservoir model

We do investigate the effects of reservoir heterogeneities through a numerical simulations and it requires a specific spatial distribution of the reservoir properties. Since geological constituent variations have spread to several scales and available data are so limited, a geostatistical approach is needed for building a geological model. In addition, stochastic probabilistic modeling is required for sensitivity purposes. Unfortunately, this approach is costly and time-consuming. As a result, there exists a great incentive to create a simpler physical model for faster prediction of reservoir performance. To do this, the percolation theory has been successfully applied in the modeling of geometrical and physical properties such as connectivity. The advantage is that it reduces many results to some master curves from which all possible outcomes can be predicted by simple algebraic transformations. These transformations mean that the effects of the complex geometry, which influence the flow, can easily be estimated, in a fraction of a second on a spreadsheet. It is emphasized that accuracy is not lost or little is lost. In particular, the percolation approach can estimate the uncertainty in these properties which is usually not possible with a Single realization model.

The disadvantage is that some of the flow physics and heterogeneity detail are missed [33], But considering the significant edges of this simulation compared to its downsides it must be considered a fast and low-cost alternative.

The model is generated by randomly distributing the sand bodies with specific geometrical shapes (such as squared sand bodies in 2D or cubes in 3D) within an impermeable background (e.g., shale). Then the sand occupation fraction is the reservoir Net-To-Gross, NTG [26]. The percolation theory was first conceptually developed by Broadbent and Hammersley in 1957 [27, 34]. But it was de Gennes’ in 1976 who first studied diffusion analysis on percolation structures by using random walks [35].

Next, it was mathematically defined using geometrical and probabilistic concepts (Stauffer and Aharony [36]). Since then this topic has been developed in many disciplines (Sahimi [37]). For example, this theory has been applied to study of the flow in porous and fractured media [37]. In particular, it proposes a framework for diffusion and dispersion of flow in porous media (Hunt et al., 2014) [38].

There are various modes of percolation problems. The simplest form is site percolation where the sites of a lattice are randomly occupied to a fraction p called occupation fraction [39]. It is assumed that the occupation of one site is independent of the occupation of any other site. The occupied sites represent the permeable part of the system and consequently the flow passes only through the permeable connected zones. Alternatively, one may occupy the bonds instead of the sites, known as bond percolation (Stauffer and Aharony, 1992 [36]). This model can be useful for modeling flow through fracture networks. Percolation structures could be made off-grid by using geometrical permeable objects that are distributed randomly in an impermeable space. This is called continuum percolation (e.g., Berkowitz [40]; King [41]).

Also, there are other percolation modes such as correlated percolation, directed percolation or dynamic percolation. In correlated percolation the objects follow some sorts of spatial correlations (Lee and Torquato [42]; Prakash et al. [43]), whereas in dynamic percolation mode the sites/bonds are not fixed and there is a probability that a bond or a site is open or closed depending on the time. Furthermore, one may set a direction (for example from bottom toward the top of the system) or the so-called directed percolation, which is appropriate for some cases of two-phase displacement (Wilkinson and Willemsen [44]).

In classical percolation, clusters of occupied sites (or permeable sites) are the main entities that control the local conductivity. As occupation fraction p is increased, clusters of permeable sites grow in size and at a specific occupancy (called percolation threshold), one large cluster becomes the controlling cluster in the medium. This critical value in 2D and for site percolation is about p_c = 0.59275 [34]. This main cluster is called spanning, infinite, or percolation cluster. If the occupation fraction p exceeds the threshold value p_c, the smaller cluster gets absorbed in the spanning cluster and create a cluster of compartments that connect the two sides of a system [26].

The other percolation property is percolation probability P (p) (connectivity or connected fraction), which is defined as the probability of any occupied site that belongs to the spanning cluster for a given value of p [33]. Below the percolation threshold, there is no connectivity while above the threshold, the percolation probability P(p) increases rapidly by the scaling law [45]. Equation (4) shows for the infinite size systems and in the vicinity of the percolation threshold (i.e., just above the threshold), a simple power-law applies [45]:

(4)

4 Results and discussion

4.1 Method of permeability map generation

To generate the permeability maps we used the MATLAB software, which is a matrix based programming medium. As the reservoir model for this study consists of 101 × 101 cells in two dimensions, we generated (using the rand command) 10 201 random numbers for the permeability property was taken from a given distribution, that 10 201 cells between 0 and 1 by random selected.

Next, the random numbers were regenerated with a normal logarithm distribution and the permeability maps were first generated with same command for 101 × 101 grid cells. Then the average permeability (in logarithmic scale), range, covariance model, Gaussian variogram, sill, standard deviation (logarithm), and variance (logarithm) were determined.

The simulation is based on the standard finite differences scheme, semi-implicit discretization, while harmonic mean of transmissivities is applied between the blocks. The presence of a strong bias associated with standard finite-difference scheme of the flow equations along the assumption of harmonic means for internodal transmissivities, results in a very slow convergence [46].

Moreover, the mean and variance of effective hydraulic conductivity distribution in the medium depend on the scale of coarse grids. This is true especially for binary media with occupation probability around the threshold [47]. To overcome this effect in our simulation, the refinement of grids with a ratio order of 10 is been used. Also, a sensitivity analysis has been performed to watch the results thoroughly.

The heterogeneous reservoirs are simulated with of blocks of two different permeabilities, 0 and a random variable k. Herein, the permeability of inactive blocks is set to be zero, and the permeability number of the block is considered if it is active. So the active blocks are chosen randomly in the model with the p proportion. Grid blocks of zero act as a barrier, neither fluid enters nor fluid leaves. In the random selection of grid blocks, the well block grid of the well must be actively selected.

First, we generate 100 independent permeability maps (realizations). Thorough these maps, the cells are randomly occupied by the fractions of p = 50%, 60%, 75%, 90%, and 100%. Good to mention that in the lower occupation fraction, for example when p = 60%, small size clusters may be formed. However, as the occupation fraction increases, the cluster size of permeable cells turns to be larger and larger. Now, the spanning cluster is the one that connects both sides of the system. For the cases of this study (site percolation), this occurs at an occupation probability around p = 59%. Hence, we generate the permeability maps (realizations of the permeability models) by considering the occupancy probabilities greater than p = 59%. The grid which hosts the wellbore (in the middle of the reservoir) is by default assumed to be occupied (with permeability k). To illustrate this matter we produced some realizations of permeability models at occupation probabilities of p = 50%, 59%, 60%, 75%, 90%, 99% and 100% by the MATLAB software are shown in Figure 2.

Fig. 2 Illustration of clusters formed at various occupancy probabilities p =  50%–100% where yellow shows the permeable cells.

4.2 Simulation details

For the purpose of numerical simulation, a 2D reservoir (100 ft, 100 ft, 200 ft) with a well located in the middle is been simulated. To compute the system properties such as occupation fraction, p, we cover the model with a fine grid that consists of 101 cells in the X direction, 101 cells in the Y direction, and one cell in the vertical direction consisting of 10 201 blocks. Therein, the production wellbore is located in the (51, 51, 1) block.

Then, these permeability maps are inserted to the flow simulator input data file, and the well constraints are allocated. Running the flow simulator for each realization of the model, the results for pressure versus time were collected. The simulation was repeated for 100 realizations at the same occupancy probability, p, to get reasonable and repeatable pressure-time data. The analysis of this pressure profile can give the effective permeability of the models.

Due to the application of flow simulations to take the well testing analysis, we need some details of the reservoir rock and fluid properties (e.g., viscosity and density for oil 0.51 cp, 37.04 lbm/ft³ and for water 0.31 cp, 62.96 lbm/ft³) which was taken from Odeh work [48] as well as reservoir thickness (200 ft), initial pressure and temperature (4800 psi, 200 °F) and well information (well radius r_w = 0.25 ft).

As the oil is assumed be a dead oil, water is the single-phase fluid with pressure-dependent properties. But, we have assumed water to be immovable in this study. On the other side and in the dead oil model, it is assumed that the density of fluids is only a function of pressure; otherwise, the default values of the simulator is used to calculate the equilibrium equation. After creating the data file in the simulator, we performed the simulations and then we extracted the time and pressure data in separately.

Well testing Software (Ecrine v4.02 software) is utilized to this simulation. To construct the well model we made the following assumptions: Reservoir temperature is 200 °F, well radius is 0.25 ft, reservoir thickness is 200 ft, porosity is 0.3, fluid is single-phase, and well storage is neglected. The results of all realizations at various occupation probabilities are shown in terms of the pressure and pressure derivative against time (Figs. 3–4). Moreover, the averaged curves for the pressure and pressure derivative against time are demonstrated in Figures 3–4.

Fig. 3 The plot of results of all realizations pressure and pressure derivative against time on the Log-log scale at occupation fraction (p = 0.5, 0.6, 0.75, 0.9, 1).

Fig. 4 The plot of average curves pressure and pressure derivative against time on the Log-log scale at occupation fraction (p = 0.5, 0.6, 0.75, 0.9, 1).

The analysis of pressure-time for each case provides the corresponding effective permeability. It is been noticed that averaging the pressure and pressure derivative curves, that yields some average of 1/K_efff, does not coincide with the average of K_eff unless the variance is very small.

The results of the analysis for pressures-times are summarized in Table 1.

Comparison of the average curves of the pressures and pressure derivative versus time on Log-log scale of these five models (p = 50%, 60%, 75%, 90%, 100%) are shown in Figure 5.

Fig. 5 Illustration of the average curves for the pressure and pressure derivative against time on Log-log scale at various occupation fractions.

Moreover, the average pressure (P) and the average pressure derivatives (P′) versus time obtained from all realizations in these five different occupation models are depicted in Figure 6.

Fig. 6 (a) The average pressure (P) and (b) the average pressure derivatives (P′) against time obtained from all realizations of five models (p = 0.5, 0.6, 0.75, 0.9, 1).

Additionally, standard deviation of pressure (Δp) and standard deviation of pressure derivatives (Δp′) versus time obtained from all realizations in these five different occupation models could be seen in Figure 7.

Fig. 7 (a) The standard deviation for pressure (Δp) and (b) the Standard deviation for pressure derivatives (∆p′) against time obtained from all realizations of five models (p = 0.5, 0.6, 0.75, 0.9, 1).

The average permeabilities calculated from these five models (p = 0.5, 0.6, 0.75, 0.9, 1) are graphed in Figure 8. As can be seen, the dependency of the average permeability to the occupancy probability follows a power law behavior. The power law exponent of 1.2 is in agreement with the critical exponent of 2D systems from the percolation theory.

Fig. 8 Illustration of the variation of average permeability (over all realizations) to the occupancy probability above the threshold which shows a power law behavior with exponent μ = 1.2.

The analysis of pressure-time of each regenerated permeability maps with the normal logarithm of permeability values, are performed to obtain the corresponding effective permeability values. The results are summarized in Table 2.

The average permeabilities of these five models (p = 0.5, 0.6, 0.75, 0.9, 1) are shown in Figure 9. Obviously, the average permeability related to the occupancy probability with a power law behavior as was the case before.

Fig. 9 Illustration of the variation of Logarithm average permeability (over all realizations) to the occupancy probability which shows a power law behavior.

5 Results and discussions

The nature of fluid flow in hydrocarbon reservoirs is considered a complex subject because of the complicated geometries resulting from complex sedimentary processes. To determine the uncertainty in reservoir performance or for sensitivity purposes, constructing several possible reservoir stochastic realizations is been a standard solution, but it is computationally time-consuming. In this research percolation is applied as a mathematical theory of connectivity which is capable of investigating the impact of the geological uncertainties on the reservoir’s performance.

Herein, we showed the geometrical and physical properties (e.g. connectivity and conductivity) could be estimated directly by using percolation theory. We generated a random binary field (i.e. the rock is either permeable or impermeable) for the permeability distribution and assumed that the flow path is strongly controlled by the connectivity of permeable units and not strongly modified by the flow dynamics themselves.

It is been found out that someone can have reliable estimates of effective permeability for very heterogeneous models, predict uncertaintly in performance parameters at field scale and to determine some properties of rock at pore scale very quickly, if they apply the percolation theory, correctly.

It was seen that an increase in the occupancy probability (i.e. the fraction of permeable cells in the permeability map) improves the effective permeability of the medium by enhancing the connectivity of permeable regions (clusters in percolation terminology). This results in a decrease in the pressure drop across the medium and consequently affect the reservoir productivity that is shown in Figures 3–8. Also, it is been demonstrated that by increasing probability p, the trend of pressure graphs becomes more accurate and sharper that could be seen in Figures 3–6. The same is true of error or standard deviation, which is illustrated in Figure 7.

Moreover, the dependency of the effective permeability to the occupancy probability follows a power law relation with a measured power-law exponent in line with the percolation critical exponents. Being so, the results of this study is consistent with the results of previous studies and therefore the modeling and computations are partly confirmed.

Finally, it is was found out that a power averaging method along using a moment from the permeability distribution (such as the geometric mean) and the permeability threshold (i.e., showing the structure and topology) could give us estimates for the effective permeability.

The result of this research opens insights on new methods in estimating the physical properties of real heterogeneous porous media by using concepts from the percolation theory. This approach can be extended to real 3D models and to consider.

The approach presented in the work is general although the simulation model used in this work and the type of well test was simple. We may perform this analysis on more complex binary porous mediums with small/large correlation ranges or strong anisotropy. Moreover, the other types of well tests may be used.

Acknowledgments

The referees are gratefully acknowledged for their useful comments.

References

All Tables

All Figures

	Fig. 1 Diagnostic plot of reservoirs. (a) Homogenous reservoir with constant boundary pressure, (b) homogenous reservoir with no-flow boundary, (c) dual porosity of the naturally fractured reservoir, (d) triple porosity of the naturally fractured reservoir [8].
In the text

	Fig. 2 Illustration of clusters formed at various occupancy probabilities p =  50%–100% where yellow shows the permeable cells.
In the text

	Fig. 3 The plot of results of all realizations pressure and pressure derivative against time on the Log-log scale at occupation fraction (p = 0.5, 0.6, 0.75, 0.9, 1).
In the text

	Fig. 4 The plot of average curves pressure and pressure derivative against time on the Log-log scale at occupation fraction (p = 0.5, 0.6, 0.75, 0.9, 1).
In the text

	Fig. 5 Illustration of the average curves for the pressure and pressure derivative against time on Log-log scale at various occupation fractions.
In the text

	Fig. 6 (a) The average pressure (P) and (b) the average pressure derivatives (P′) against time obtained from all realizations of five models (p = 0.5, 0.6, 0.75, 0.9, 1).
In the text

	Fig. 7 (a) The standard deviation for pressure (Δp) and (b) the Standard deviation for pressure derivatives (∆p′) against time obtained from all realizations of five models (p = 0.5, 0.6, 0.75, 0.9, 1).
In the text

	Fig. 8 Illustration of the variation of average permeability (over all realizations) to the occupancy probability above the threshold which shows a power law behavior with exponent μ = 1.2.
In the text

	Fig. 9 Illustration of the variation of Logarithm average permeability (over all realizations) to the occupancy probability which shows a power law behavior.
In the text