Million node fracture: size matters?

Transmissivity of self-affine fractures was computed numerically as a function of the grid size. One-million-node fractures (1024 × 1024 nodes) with fractal dimensions of 2.2–2.6 were cut into successively smaller fractures (“generations”), and transmissivities computed. The number of fractures in each generation was increased by a factor of 4. Considerable scatter in transmissivity was observed for smaller grid sizes. Average transmissivity of the fractures in the generation decreased with the grid size, without approaching any asymptotic value, which indicates no representative elementary volume (REV). This happened despite the average mean aperture being the same in each generation. The results indicate that it is not possible to estimate the transmissivity of a large fracture by cutting it into smaller fractures, running flow simulations on those and averaging the results. The decrease in transmissivity with the grid size was found to be due to an increase in the flow tortuosity.


Introduction
Fractures play an important role in controlling hydraulic and transport properties of fractured rocks. In some types of hydrocarbon reservoirs, they can largely control the reservoir flow (Aguilera 1980;Li et al. 2021). Understanding flow, transport and fluid displacement in natural fractures is crucial for optimization of production from such reservoirs. Understanding fracture flow and fracture permeability is also important for optimization of drilling operations and drilling fluid design since fracture properties often control mud losses (Majidi et al. 2010). Fractures also control groundwater flow and contaminant transport in certain types of aquifers (Maldaner et al. 2021).
Experimental studies of fluid flow in fractures should, ideally, be conducted in situ. However, due to substantial practical difficulties associated with such experiments, permeability of rough-walled fractures has often been studied analytically (Lenci et al. 2020;Rodríguez de Castro and Radilla 2017;Wrobel et al. 2021) and numerically (Zhang and Huang 2018), using synthetic, numerically generated fractures or using digital replicas of real fractures. Fracture surfaces are known to have self-affine landscapes, i.e. the surface looks essentially the same if we zoom in or out on it. To describe the geometrical properties of such surfaces, the fractal dimension, D, is used (Mourzenko et al. 2018). Fractal dimension is a measure of surface roughness or "heterogeneity"; higher fractal dimension results in greater flow tortuosity. Other parameters describing the aperture map of a fracture are its mean aperture, w m , and the root mean square (rms). The mean aperture is simply the arithmetic average of all local aperture values. The rms is the standard deviation of the local aperture from w m .
Numerical studies of fracture flow have been performed since at least the work of S. R. Brown (Brown 1987). In a typical numerical simulation of this kind (Fig. 1), a rectangular or square fracture plane with a distribution of aperture is discretized using a regular Cartesian grid in the fracture plane (x, y). Lubrication theory approximation, represented by the Reynolds equation, is then used to approximate the average fluid velocity locally (Zimmerman et al. 1991). The pressure field in the fracture is then obtained by solving the continuity equation numerically, e.g. using the finite-difference method. In Eq. (1), L is the size of the square fracture; w is the local fracture aperture; P is the fluid pressure. We assume cubic law in Eq. (1) since this is by far the most popular closure law used in fracture flow modelling today, even though its validity has been duly challenged several times (Cvetkovic et al. 2004;Mourzenko et al. 1995;Xu et al. 2021). The boundary conditions are of Dirichlet type on the left and right edges of the fracture, and of Neumann type (no flux) on the top and bottom edges. The overall flow direction is thus along the x-axes. Once the solution has been obtained, the hydraulic aperture of the fracture can be evaluated. The hydraulic aperture of a rough-walled fracture is the aperture of a fracture with smooth parallel walls that produces the same flow rate under the given pressure difference, (P 1 -P 2 ), as the rough-walled fracture. Hydraulic aperture, w h , can be either smaller or larger than the mean aperture of the fracture (Méheust and Schmittbuhl 2000). Fracture permeability is equal to w 2 h 12 . Fracture transmissivity is defined as w 3 h here. Numerical studies of fracture flow have been mostly performed using grids of 10 3 -10 5 nodes in the plane of the fracture (Table 1). It is not clear how representative such studies are for hydraulic behaviour of real fractures. So far, few studies have addressed the size effect on hydraulic properties of rock fractures. In (Raven and Gale 1985), specimens of different sizes were made from different parts of the same natural fracture. Fracture transmissivity was found to decrease with increasing sample size. This was attributed by the authors to greater flow tortuosity in larger fractures. In (Koyama et al. 2006), no clear trend in the fracture transmissivity was observed when the fracture size was increased from 50 × 50 mm to 100 × 100 mm and to 200 × 200 mm, unless the fracture faces were sheared relative to each other. In another numerical study (Matsuki et al. 2006), a perhaps expected conclusion was reached that when a characteristic scale (maximum spatial wavelength) is introduced to the fracture aperture distribution, the size effect disappears at fracture sizes above that cut-off value. The cut-off size is the fracture size above which the rms of the aperture remains constant (instead of increasing with fracture size, as would be the case for a self-affine aperture map). Such characteristic length may be found in fracture landscapes built using parameters obtained from artificial fractures created in laboratory tests, where the fracture surfaces are well matched at larger spatial wavelengths, as was the case in (Matsuki et al. 2006). Another assumption made in (Matsuki et al. 2006) was single-point contact between fracture surfaces. A characteristic length of the aperture distribution was discovered in some natural joints when their faces were artificially positioned so as to match each other (Brown et al. 1986). In other cases, however, especially with mineralized or weathered surfaces, and with nonzero shear displacement between them, the "matedness" of the faces might be vanishing even at large scales.
If a representative volume does exist, the hydraulic aperture should be converging to a certain value as the grid size increases. The objective of our study is to investigate whether fracture transmissivity converges to a certain value with the grid size. Another issue we are interested in is whether small-grid simulations could be used to predict the transmissivity of a large, million-node fracture. For instance, is transmissivity of a 1024 × 1024 fracture approximately equal to the average of hydraulic apertures of four 512 × 512 fractures cut out of this 1024 × 1024 fracture? If the answer is positive, then performing these four simulations would provide a much faster way to determine the hydraulic aperture of the 1024 × 1024 fracture.

Simulation setup
The key idea was to generate a synthetic fracture of 1024 × 1024 nodes ("one million node fracture"), compute its transmissivity, then successively subdivide this fracture into smaller fractures, compute their transmissivities, etc. The million-node fracture was generation zero (g = 0). After computing its hydraulic aperture, the fracture was cut into 4 equal fractures (generation 1, g = 1). Their hydraulic apertures were computed, and then this process, shown schematically in Fig. 2, was recursively repeated down to fractures of 32 × 32 nodes (g = 5). The number of simulations in each generation g was 4 g . For each generation, average hydraulic aperture was calculated, along with some other parameters described below. Five million-node fractures were generated using the spectral synthesis method (Huang et al. 1992). The fractal dimension of this initial aperture map was verified by the variogram technique (Ruello et al. 2006) and by an improved version of the box-counting technique described in (Ai et al. 2014). The fractal dimensions, D, set in the spectral generator were 2.2, 2.3, 2.4, 2.5 and 2.6 for aperture maps A, B, C, D and E, respectively. These are within the range usually quoted for rock fractures (Kulatilake and Um 1997;Lanaro 2000;Odling 1994;Poon et al. 1992;Power and Durham 1997;Schmittbuhl et al. 1993). The fractal dimensions obtained from the variogram analysis were found always to be equal to the dimensions set in the generator. The boxcounting dimension was typically slightly smaller, by up to 0.1. These five initial aperture maps had mean zero aperture. In order to create the zero-generation fractures from the initial aperture maps, the maps were rescaled by the same factor so that the percentage of overlaps (contact spots) be within 3 to 9% when 1 mm is added to the initial map. The percentage of overlaps was higher in fractures with greater fractal dimension. The five rescaled aperture maps were then used to generate the five million-node fractures, by adding 1 mm to the rescaled map. In the resulting aperture map, nodal apertures smaller than 0.1 mm were set equal to 0.1 mm in order to ensure convergence of the flow solver. The percentage of nodes with aperture set to 0.1 mm was from 4% (fracture A, D = 2.2) to 12% (fracture E, D = 2.6).
Two procedures were used to construct the aperture maps for fracture generations g > 0 (Fig. 3). (One million node fractures were generated in the same way in both procedures I and II.) In procedure I, the fractures of g > 0 were simply successively cut from the actual aperture map of the 1024 × 1024 fracture. In procedure II, fractures were cut from the rescaled 1024 × 1024 aperture map of zero mean aperture. Each individual map of g > 0 was offset so as to make its mean aperture equal to zero. Then, 1 mm was added to the aperture map of each g > 0 fracture, and apertures smaller than 0.1 mm were replaced with 0.1 mm. Procedure I was applied to fractures A to E. Procedure II was applied only to fracture A and was found to produce results virtually identical to procedure I (see Sect. 3). Therefore, for other fractures, only procedure I was used. An example of aperture map for a generation 3 fracture (128 × 128 nodes) obtained from the original one-million-node fracture of D = 2.4 using procedure I is shown in Fig. 4.
Steady flow simulations were performed for all fractures using finite-difference method with asymmetric scheme and arithmetic averaging of transmissivity at cell faces (van Es et al. 2014). Stabilized biconjugate solver was used; the solver tolerance was set to the same value in all simulations. This value was chosen so as to ensure that the mass balance error (see below) be less than 1% in all simulations. In total, 8189 flow simulations were performed (five fractal dimensions × 4 g fractures in each generation with procedure I; one fractal dimension × 4 g fractures in each generation with procedure II). In each simulation, the following data were calculated: 1. Mean aperture, w m , calculated as the arithmetic average of nodal aperture values. 2. Hydraulic aperture, w h . First, average discharges through the right (Q r ) and left (Q l ) boundaries of the fracture were computed: Q = (Q r + Q l )/2. Then, w h was obtained as the aperture of a smooth-walled fracture that, given the same pressure difference, would produce a discharge equal to Q. 3. Mass balance error, defined as |Q r -Q l |/Q. 4. Flow tortuosity, τ, evaluated as the ratio between the mean value of the fluid velocity and the mean value of the x-component of the fluid velocity (Koponen et al. 1996). It was shown earlier that τ is capable of capturing e.g. channelling effects in fractures (Bao et al. 2017). From the above data, the following parameters were evaluated for each generation and each fracture construction procedure: 1. Average hydraulic aperture for fractures, < w h > . 2. Ratio of the average values of hydraulic and mean apertures, < w h > / < w m > .
Averaging was done over the 4 g fractures in each generation.

Results
Hydraulic aperture averaged over 4 g simulations in each generation, g, was found to decrease with the grid size (Fig. 5). There was no tendency to approach a representative elementary volume (REV) limit, at any fractal dimension. Average hydraulic apertures obtained on smaller fractures could not be used to predict the hydraulic aperture of the million-node fracture. This outcome was obtained with both procedure I and II (Fig. 6).
Despite the clear descending trend evident in Fig. 5, there is a significant scatter in hydraulic aperture values obtained in individual simulations in the same generation, see Fig. 7 where we plot results of all flow simulations obtained with fracture C.
The decrease in w h with the grid size in Fig. 5 takes place despite the fact that the average value of the mean aperture Rescale the initial aperture map (all by the same factor) so that, if 1 mm is added to the map, the percentage of overlaps be within 3-9 %.
Add 1 mm to the rescaled map.
Set apertures <0.1 mm equal to 0.1 mm.
Cut g = 0 fracture into g = 1 fractures Procedure I for g > 0 Cut g = 0 fracture into smaller fractures Offset to make w m = 0 for each smaller fracture.
Add 1 mm to each smaller fracture.
Set apertures <0.1 mm equal to 0.1 mm.

Procedure II for g > 0
Initial map

Rescaled map
One million node fracture Cut g = 1 fracture into g = 2 fractures Continue until g = 5 in each generation, g > 0, is the same and is equal to that of the million-node fracture, by construction (cf. details of procedure I in Sect. 2). In order to find the reason for the decrease in < w h > with the grid size, it is instructive to look at the average flow tortuosity vs. grid size. Figure 8 shows that the flow tortuosity increases with the grid size, and this increase is more pronounced for fractures with higher D, as expected. The increased tortuosity is responsible for the descend in < w h > in Fig. 5.

Discussion
Our findings provide an answer to the question asked in the title of this article: size does matter, and even one million nodes are not enough. It means that most of the numerical studies performed so far on synthetic (i.e. numerically generated) fractures or numerical replicas of real fractures are, strictly speaking, not representative. This conclusion does not invalidate their results (or those of future studies), nor does it make them less valuable. It only shows that care should be exercised in their interpretation and, particularly, in extrapolating these results to the in situ scale. We, therefore, suggest that the choice of the grid size in fracture flow simulations be based on considerations of numerical accuracy (e.g. the ability of the model to resolve a specific physical effect, e.g. fingering) rather than by the modeller's ambition to perform a simulation representative of the fracture flow in situ. As pointed out in (Méheust and Schmittbuhl 2000), "no small isotropic elementary volume can be defined for the hydraulic properties". Our results suggest that no elementary volume can be defined for hydraulic properties that could be a standard (or even a guideline recommendation) for numerical and laboratory studies of flow and transport in rock fractures. A million-node fracture is currently a practical limit that can be run on a desktop computer. Our results suggest that this is by far not enough to call such simulations representative. In (Brown et al. 1986;Matsuki et al. 2006), decrease in transmissivity with the fracture size was attributed to an increase in the aperture rms (standard deviation). Our results suggest that this might not always be the case: In Fig. 5, the hydraulic aperture keeps dropping even after the rms begins to stabilize. The reason for this trend in w h is increase in tortuosity that continues as the fracture size increases. Our results confirm the conjecture made in (Raven and Gale 1985) that larger tortuosity of flow in larger fractures is responsible for their reduced hydraulic aperture (Fig. 8).
In (Raven and Gale 1985), it was indirectly supported by laboratory measurements of flow and deformation. It was speculated that in smaller fractures, there are only a few high asperities; in larger fractures, there are more asperities, and they are more evenly distributed over the fracture plane. Thus, the tortuosity in larger fractures must be greater and the transmissivity smaller.
An extra series of simulations was set up with a Gaussian fracture, i.e. with the aperture at each node simply sampled from a Gaussian distribution. The average aperture of the million-node fracture was 1 mm, rms 0.7 mm, percentage of contacts 10%. The averaged hydraulic aperture, < w h > , remained virtually constant and independent of the grid size in this case. Thus, in order to evaluate the hydraulic aperture of the million-node fracture in this case, the fracture could be cut into e.g. 32 × 32 square fractures, flow simulations could be run on each of them, and the average w h could be used as an estimate for the million-node fracture. Interestingly, the average flow tortuosity increased from 1.0519 to 1.0535 as the grid size increased from 32 × 32 to 1024 × 1024. This increase was, however, insignificant compared to that observed with self-affine fractures (cf. Fig. 8).
Our results, in particular the decrease in the transmissivity with the grid size, are consistent with (Mourzenko et al. 2001) albeit a different procedure for constructing the subsequent fracture generations was used there. On the other hand, in discrete fracture network (DFN) models, a positive correlation between fracture size and transmissivity of individual fractures is often assumed, see e.g. (Hyman et al. 2016) and references therein. Our results suggest that on hydrodynamics grounds, transmissivity should either decrease with fracture size (for self-affine fractures) or remain constant (for Gaussian fractures). The type of correlation (positive, negative or none) affects the overall permeability of the fracture network and greatly influences the breakthrough time (Hyman et al. 2016). Further work is needed in order to understand how the assumed positive correlation between transmissivity and fracture size could be made consistent with the principles of hydrodynamics. It should be reminded here that average mean aperture in our simulations, < w m > , remained constant at all g.
Considerable scatter in hydraulic aperture values between the fractures belonging to the same generation (Fig. 7) might be responsible for the apparent lack of trend in w h vs. fracture size reported in (Koyama et al. 2006). The scatter of transmissivities at smaller fracture sizes in Fig. 7 is so large that unless the number of specimens is sufficiently large, it is unlikely that there will be discovered any trend if only a few fractures of each size are tested.
Our failure to discover REV is hardly surprising, given self-affinity of rock fractures (Schmittbuhl et al. 2008). Our results show that the averaged hydraulic aperture keeps changing as the size of the fracture increases towards a million nodes. Averaging the hydraulic properties obtained on a large number of smaller fractures cannot give an estimate of the million-node fracture's hydraulic aperture.
In our study, the fractal dimension was assumed to be independent of the fracture size. There is some evidence suggesting that the fractal dimension, in reality, may increase with the fracture size (Koyama et al. 2006). This effect might even further reduce the transmissivity of larger fractures.

Conclusions
Numerical simulations show that hydraulic aperture of rough-walled, self-affine fractures keeps decreasing almost linearly as the grid size (number of nodes along the fracture side) increases. Even a one-million node fracture (1024 × 1024 nodes) does not provide a representative elementary volume for hydraulic aperture. This is as expected for self-affine fractures. This decrease in hydraulic aperture is due to increasing flow tortuosity and not due to increasing aperture rms since the latter tends to level off as the grid size increases in our simulations.
Funding No funding was received.

Conflict of interest
There are no conflicts of interest and no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will 1 3 need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.