Abstract

Digital rock twins are widely used to obtain hydraulic properties of porous media by simulating pore-scale fluid flow. Multifractal characteristics of pore geometry and flow velocity distribution have been discovered with two-dimensional (2D) images and three-dimensional (3D) models, whereas the dependency of results on the resolution is not well known. We investigated resolution-dependent multifractal properties of 3D twin models for a sandstone sample with originally 3 μm resolution images. 3D pore-scale water flow was simulated with the lattice Boltzmann method (LBM). As indicated by multifractal analyses, the generalized dimension spectra, the Hölder exponent spectra, and singularity spectra of the flow velocity are similar to that of the pore geometry but different in ranges and sensitivities to the change in the model resolution. Nonlinear dependencies of 2D/3D porosity, holistic/slice permeability, equivalent pore radius squared, and multifractal parameters on the resolution were discussed.

1. Introduction

With the development of scanning and imaging techniques, a large number of digital rock twins have been obtained in the world to represent the pore space of rock samples. It provided abundant opportunities to reconstruct and study their hydraulic properties, or even other physical, chemical, and biological properties [1]. In general, a digital twin of a rock sample is composed of pixels in two-dimensional (2D) models or voxels in three-dimensional (3D) models signed as pore or solid. The seepage behaviours in porous media can then be simulated as fluid flow at the pore scale with the computational fluid dynamics (CFD) approach, solving the Navier-Stokes (NS) equation directly [2, 3] or applying the lattice Boltzmann method (LBM) [410]. The simulation results of the flow field in the digital twin of a rock sample are used for estimating the permeability or further modelling on solutes transport.

It has been well known that most of the porous media have fractal geometry [1113]. Multifractal analysis was widely used to reveal the fractal properties on the pore size distribution [14, 15] or the 2D/3D pore structure [1624]. Both box-counting method and the sandbox method were implemented in multifractal analyses of the pore structure [25]. A natural consequence of the fractal pore space is that the pore-scale flow velocity may also have a fractal distribution. This has been investigated from LBM results in fractal pore models [2628] that synthetically constructed from realization of a 3D random fractal lattice [29]. The multifractal distribution of the flow velocity was revealed, but the results depended on the counting method. It seems that the sandbox method is more reliable than the box-counting method [28]. However, the fractal behaviours of the flow velocity in real rock samples have not been investigated.

For digital rock twins, a few studies investigated the resolution effect on estimated multifractal spectra in the pore space. Jouini et al. [30] found that the multifractal dimension and porosity of carbonate samples varied with the magnification effect on SEM images, but the correlation was unclear. Karimpouli and Tahmasebi [20] conducted multifractal analyses on the pore geometry of sandstones and carbonates using multiresolution digital twins, and they suggested that multifractal results should be interpreted based on image resolution. However, they did not find a significant scale effect. Shah et al. [31] checked the effect of the image resolution for carbonate rocks and found that the porosity and permeability may decrease with coarsening resolution. This effect may also exist in digital twins of sandstone samples, as reported by Guan et al. [32]. In an investigation of Saxena et al. [33], the permeability estimated from coarse resolution was higher than that estimated from fine resolution. Nevertheless, the resolution effect on the simulated pore-scale flow was still poorly known. The resolution-dependent permeability of the multifractal pore structure may be linked with the multifractal flow velocity distribution, whereas the relationship is unclear.

In this study, we investigated the resolution-dependent multifractal features of both the pore geometry and flow velocity in 3D digital models of a sandstone sample, aiming to check different dependencies and potential relationship between the multifractal pore geometry and flow velocity. Nine twin models were constructed with different magnification levels on the basis of an original high-resolution image set. The simulation of the steady-state water flow at the pore-scale was implemented with LBM, assuming a fluid pressure difference on sides. The permeability of these twin models was then identified from Darcy’s law and compared to the estimated porosity. Multifractal analysis of the 3D mass distribution was performed using the sandbox method to obtain the multifractal spectra of pores and flow velocities. These multifractal spectra were observed and compared with that in existing studies to discuss the different resolution-dependent behaviours between the pore and velocity spaces.

2. Materials and Methods

2.1. Materials

In this study, a sandstone sample with a microstructure containing both intergrain and intragrain pore spaces is investigated. Digital image data of the sample [34, 35] were provided on the Digital Rocks Portal, including 3D images of slices with a voxel length of 3 μm. A cubic box with half the size of the imaged sample was extracted by cropping a square region of interest (ROI) covering pixels of the original pixels in each slice (Figure 1) among the first 256 slices. Thus, there are voxels in the cubic sample volume characterized by a physical length of 768 μm.

Digital twin models of different resolutions were built for the studied cubic sample by a magnification approach with the bicubic interpolation in ImageJ [36]. In comparison to other interpolation methods [37], such as the nearest neighbour and bilinear interpolations, the bicubic interpolation can yield smoother resampling images [38]. This interpolation was accomplished using the Catmull-Rom spline [39] in ImageJ. There are 2563 voxels in the first model, denoted as M256, with the original resolution of 3 μm. Other twin models were obtained from M256 by magnification using the bicubic interpolation with resolutions from 4 μm to 48 μm, as listed in Table 1. The digital colour was transferred from the RGB signal to the 8 bit grey value in the bicubic interpolation, varying between 0 and 255, to facilitate the binarization process.

The identified pore structure of a porous medium depends on the colour threshold to separate pores and solids [31]. Numerous methods could be used in image segmentation [33, 4043]. We checked several methods, such as the midvalue method and Otsu’s algorithm [44]. Reasonable results are required for a suitable method: when the resolution comes to be finer, the porosity converges to that of the original digital rock. As indicated, Otsu’s algorithm soundly satisfied the convergence and was finally chosen in this study. Performance of Otsu’s algorithm has also been verified in other existing studies [33, 43, 45]. In further analyses, the M256 model remains the original binarized images, while the other models are composed of images magnified with the bicubic interpolation and binarized with Otsu’s algorithm. This binarization process yields images only containing black (255, representing pores) and white (0, representing solids) voxels, as shown in Figure 1(b) for M256 and Figure 2 for other models.

After binarization, the volumetric porosity (3D porosity) was calculated as the proportion of voxels signed as pores (black) (Table 1). With the model resolution worsening (from 3 μm for M256 to 48 μm for M16), this 3D porosity slightly fluctuates between 0.328 and 0.335 when the resolution is smaller than 16 μm and then decreases rapidly with the increase in the resolution, as shown in Figure 3(a). We also calculated the surface porosity (2D porosity) on each slice ( plane) perpendicular to the direction. This 2D porosity varies with the distance along the -axis, as shown in Figure 3(b) with the normalized distance from 0 to 1. The fluctuating patterns of the 2D porosity are similar for different twin models, approximately ±0.13 around 0.36 for M256 to M48, whereas the gap between models becomes significant for M32, M24, and M16. M256 shows a high-frequency fluctuation with a small amplitude in addition to the curve of M48. In addition, the outline was recognized for each pore on a slice, and the perimeter and enveloped area were estimated from ImageJ. Results were applied to estimate the equivalent hydraulic radius () of each individual pore as a circle. The radii squared (, Table 1) were prepared for further investigation on the relationship between the permeability and pore structure (Section 3.1.1).

2.2. Lattice Boltzmann Method
2.2.1. Lattice BGK Model

The flow at the pore scale was simulated with LBM on the basis of the lattice Bhatnagar-Gross-Krook (BGK) model [46]. In LBM, the fluid flow is regarded as a huge amount of particles following a distribution function, . The distribution function describes the percentage of particles in a certain location () at a given time () within a certain range of velocities. The lattice BGK model is a single relaxation time model that mostly adopted in LBM to simulate streaming and collision of particles.

The BGK model estimates the change in the distribution function versus the lattice spacing, , and the time step, , from the following equation [5]: where is the dimensionless relaxation time towards the equilibrium distribution, ; is the frequency of the relaxation coefficient, determined as ; and the subscript character denotes the discrete direction of the velocity.

The equilibrium distribution function is determined by where is the weighting factor, is the fluid density, is the fluid flow velocity, is the unit vector along the flow direction, and is the sound speed. When , the value of is [28, 47].

Different linking schemes are developed for the lattice BGK model to arrange the connection between neighbouring nodes, denoted as , in which is the number of dimensions and represents the number of particles considered [47]. For 3D flow simulation, D3Q15 and D3Q19 are widely used schemes. In this study, D3Q19 was chosen, in which 19 velocity vectors exist, including 18 moving directions and one still (Figure 4). The corresponding weighting factors, , are 1/3 for , 1/18 for to , and 1/36 for to [28].

The fluid density, , is determined through the mass conservation of particles for each lattice: Note that is counted from 0 to 18 for the D3Q19 model. Similarly, the flow velocity vector, , can be calculated based on the momentum conservation as follows:

The kinematic viscosity, , is determined by

2.2.2. LBM Simulation on Digital Rock Twin

The pore-scale flow in a digital twin model was simulated with the LBM, covering all voxels in the model. Each voxel in M256 was arranged as a node in the lattice BGK model so that the numerical resolution of LBM simulation is 3 μm. This lattice grid was also used for other twin models, in which a voxel covered a group of nodes with the same solid or fluid phase. Estimation was active only on the pore nodes, whereas the fluid flow variables were set to zero for all solid-phase nodes. At halfway between the solid-fluid boundary, the bounce-back scheme [48, 49] was applied, meaning the “half-way wall bounce-back” scheme [50]. This scheme has the second order accuracy and is an alternative to the “complete” bounce-back scheme (bounce-back rule applied at the boundary of solids) that only has the first order accuracy [50]. The simulation started from an initial condition, in which a zero velocity was assigned for all nodes. The flow was driven by a pressure difference (0.01 Pa) on two sides at the locations of normalized and , triggering a macrovelocity of water flow in the positive direction. Key parameters of the lattice BGK model are provided in Table 2. The LBM simulation was implemented via Palabos (parallel lattice Boltzmann solver) [51] which had been released as a C++ code library. Results of the final steady-state flow are saved for further analysis.

A VTI (Visualization Toolkit Image Data) file was exported from the LBM simulation for a model, with the velocity and pressure data of all nodes. The software, ParaView, was used to process the output results, which could display the 3D velocity fields.

According to Darcy’s law, the permeability, , can be obtained through the parameters in Table 2 and the LBM results, using the following equation: where is the Darcy velocity along the direction and is the pressure gradient. This calculation can be performed in the whole (holistic permeability) or along a twin model slice at a special distance (slice permeability). In calculating the holistic permeability, is estimated as the average value of for all nodes (including zero value at solid-phase voxels), and is 13.02083 Pa/m as given in Table 2. For the slice permeability, is estimated from the data of nodes on a model slice at , and is approximately estimated with the central difference method, i.e., estimated from the pressure difference of neighbouring slices.

2.3. Multifractal Analysis

Multifractal analysis [52] provides a singular set, , defined in a fractal structure, which is composed of subsets for probabilities with respect to different singular scale component, .

2.3.1. Sandbox Method

The sandbox method [27] considers the mass (sum of measurements, for example, the pore volume or the absolute velocity), , in each cell (a 3D cube) of a given size, . The probability mass distribution is calculated by where is the mass in cell of an initial size (smallest), , and is the number of cells in the structural model with respect to the initial cell size. The partition function, , of the order , is then estimated as where is the number of cells with respect to the size of , serves as a “microscope” in exploring different regions of the singular measurement [52], and is the mass exponent defined as

In practice, however, is approximately estimated for each by fitting the log-log plot of data points versus with an optimized slope [53].

The generalized fractal dimension, , is calculated from the mass exponent function:

When , the L’Hospital’s rule [54] is applied, which leads to

is the information dimension (or entropy dimension), quantifying the degree of measurement distribution heterogeneity. and are the capacity dimension and the correlation dimension, respectively. The plot of versus is called the generalized dimension spectra.

The Hölder exponent, , and the Hausdorff dimension, , can also be determined from the mass exponent by the Legendre transformation [55, 56]:

The plot of versus is called the multifractal spectrum or singularity spectrum.

2.3.2. Targets of Multifractal Analysis on Digital Rock Twin

The multifractal properties of different twin models were analysed for the pore geometry and flow velocity distribution.

For the pore geometry of a studied model, the mass measurement was 1 on a pore voxel, while it was 0 on a solid-phase voxel. For the velocity distribution, the mass measurement was specified as the velocity magnitude, , at a node in the lattice grid, as follows: where , , and were flow velocities along the , , and directions, respectively, obtained from the LBM simulation. For each twin model, a 3D flow velocity matrix (data of ) on all LBM lattice nodes (2563) was prepared for the multifractal analysis.

We developed a C++ program to estimate the partition function, , with Equations (7) and (8), for several specified values of . The initial cell size, , was specified as the resolution of a digital twin model for the pore geometry or the lattice spacing for the velocity distribution. The cell size, , increased step by step with an increment that was several times of , until it reached the side length of the model. We also developed a MATLAB program to obtain by fitting the log-log plot of data points versus linearly with an optimized slope. Then, the generalized dimension spectra, , and the multifractal singularity spectra, F(α) versus α(q), were estimated through Equations (10), (12), and (13), where the value of was specified from -10 to 10 with an interval of 0.2.

2.3.3. Processing Steps and Notes

The multifractal analysis in this study was implemented through three steps. First, check the data of the mass measurement to ensure that a multifractal characteristic exists. Second, estimate the multifractal spectra for a digital twin model. Third, compare multifractal parameters among the nine models of different resolutions.

To verify whether the pore structure or velocity distribution is multifractal, linear relationship between the logarithmic partition function versus logarithmic was checked; the capacity dimension, , the information dimension, , and the correlation dimension, , were compared with general experiences. A multifractal pattern exists only when , with linear relationship between and . Furthermore, the multifractal distribution usually shows a single-peak curve for the plot of versus . When almost keeps a constant value, it would not be a multifractal distribution.

Additional parameters are introduced to quantify the changes in the dimension and the Hölder exponent with when the value is larger or smaller than zero, as shown in the following:

The width of the left part () of the generalized dimension, defined by , and that of the Hölder exponent defined by , which is also the width of the right branch of , can reflect the low probability density of mass (i.e., dominance of small pore or low velocity). In comparison, the widths of the right part () of and , defined by and (also the left branch of ), respectively, correspond to the higher probability density. Therefore, multifractal parameters , , , and can be used to indicate the pore or velocity distribution in lower probability areas; meanwhile, , , , and can be used for higher probability areas.

The digital twin models of the sandstone with different resolutions can be characterized by different values of a special multifractal parameter. It was found in [12, 17] that a smaller indicates greater heterogeneity. A lower may indicate a higher degree of local mass density. Higher values of and may also indicate more complexity in the mass distribution, as well as the heterogeneity of porous media.

3. Results

3.1. LBM Simulation Results
3.1.1. Darcy Velocity, Permeability, and Equivalent Pore Radius

The results of the LBM simulation for all digital twin models were exported in VTI files, including the flow velocity and pressure data on all lattice nodes. The velocity data list , , and , as well as calculated with Equation (14). The statistics of the flow velocities are shown in Table 3. Although the macropressure gradient acted on the direction, it generated nonzero and values at most of the pore nodes. The average value is negative for all twin models. The average value is positive or negative and closer to zero than the average value. The flow velocities along , , and have similar variation levels, as indicated by the standard deviation data.

The holistic permeability was calculated according to Equation (6) with the holistic Darcy velocity, . The average velocity magnitude and the holistic Darcy velocity are shown in Figure 5(a) for comparison. As indicated, the average value is higher than the average value because of nonzero and values, while both of them almost keep constant first and then increase with the resolution coarsening, exhibiting a nonlinear response. The holistic Darcy velocity is about 1/3 of the average value of the pore nodes. The changes in the holistic permeability and 3D porosity with the increasing resolution are different, as shown in Figure 5(b). They both show slight fluctuations around the original value of M256 when the resolution is smaller than 10 μm. However, when the resolution coarsens from 16 μm (M48) to 48 μm (M16), they show opposite response: the porosity decreases while the permeability increases. This would be an unexpected phenomenon if the positive correlation between porosity and permeability was only considered. It has also been reported by Xia et al. [57] that using only the porosity could not completely explain the change in permeability for complex structures. Therefore, we need to observe another control of the permeability: the size of pores in average. This was the motivation of estimating the equivalent hydraulic radius (). The data of radii squared () in average are listed in Table 1 and plotted in Figure 5(b) in a logarithmic form of the value related to M256. It exhibits a significant increasing pattern with the coarsening process when the resolution is larger than 10 μm, i.e., more large pores were recognized with coarser images. We argue that this might be a key cause of the increase in the estimated permeability, but the effect is slightly weakened by the decrease in the 3D porosity. Quantitative analysis on the relationship will be further presented in Section 4.

We also calculated the slice permeability using Equation (6) (Section 2.2.2) and compared the results with the 2D porosity for different twin models. Figure 6 shows typical results of M256, M128, M64, and M32. The value was too small to obtain a reasonable result when was close to 0. Therefore, the permeability of some slices near the upper-stream boundary is not provided. In Figure 6, the slice permeability is exhibited with l because the data cover several orders of magnitude. For each model, the slice permeability shows a fluctuation pattern along the distance that is similar to the 2D porosity but has different peak/valley locations. The 2D porosity distribution curves for different models are similar to each other, as indicated in Figure 3(b). In comparison, the fluctuation curve of the slice permeability changes significantly with the varying resolution, as indicated by the significant difference between that of M256 (Figure 6(a)) and M32 (Figure 6(d)).

We checked the potential representative elementary volume (REV) of the studied sandstone sample. It was found that the REV size is most probably close to 800 μm, being slightly larger than the size of twin models (768 μm). Therefore, Darcy’s law may be not rigorous for the flow in a twin model. The slice and holistic permeabilities estimated herein could be regarded as the apparent or equivalent permeability for a volume that is smaller than REV.

3.1.2. 3D Distribution of the Velocity Magnitude

The velocity magnitude () fields obtained from the LBM simulations on the pore-scale flow through the nine models are shown in Figure 7. These 3D distribution patterns were displayed in ParaView using the imported VTI files. It is clear that the flow velocity has a significantly heterogeneous distribution. Preferential flow passageways exist in the space. High-velocity flow ( m/s) mainly occurs in the upper part (with relatively large value) of this sandstone sample, and the passageways mainly extend along the direction. The other parts of the flow space are tubular low-velocity zones extending vertically or horizontally. Branches of the preferential flow passageways almost follow the same pattern when the model resolution becomes coarse from M256 (Figure 7(a)) to M48 (Figure 7(f)), whereas change the shape and stretch directions when the model changes from M48 (Figure 7(f)) to M16 (Figure 7(i)). Those preferential flow passageways are produced by the heterogeneous connection of pores as exhibited in Figure 2. For a pore tube, the flow velocity is maximum in its centre, forming a circular distribution of . This pattern can be clearly identified in cross-sectional views of the velocity magnitude, as particularly shown in Figures 8(a)8(c), for M256, M48, and M16, respectively.

3.2. Multifractal Characteristics
3.2.1. Multifractal Spectra of the Pore Geometry

The log-log plots of the partition function, , when and , are typically shown in Figure 9 for the pore geometry in different twin models. A strong linear relationship between and is exhibited for decades of the cell size in sandbox counting, which indicates that this sandstone sample is a fractal porous material. Note that this linear relationship is negative for and positive for .

The graphs of the generalized fractal dimension, , the Hölder exponent, , and the singularity spectra, , are shown in Figure 10. Results of major multifractal parameters are listed in Table 4.

As shown in Figure 10(a), is negatively dependent on in a nonlinear manner. Data in Table 4 indicate that both and slightly increase with the worsening of the resolution from M256 to M48 and then decrease more from M48 to M16, while monotonically and slightly decreases with increasing resolution. They almost keep steady in the range between 2.6 and 2.9, implying a less-than-3.0 dimension in the geometry. Figure 10(b) shows nonlinear dependencies of the and values on the resolution, where the minimum value exists at 16 μm (M48 among the models). In contrast, a monotonically positive correlation between the value and the resolution is exhibited (Figure 10(b)). The variation patterns of are similar to that of , as indicated in Figures 10(c) and 10(d), while the nonlinear curve of looks more like a reverse-S shape. The Hausdorff dimension, , is a function of the Hölder exponent, , as shown in Figure 10(e). Even so, it ultimately depends on because depends on , and Figure 10(f) shows the relationship between and , which is also nonlinear.

For all of the digital twin models, the linear relationship between and is significant as that shown in Figure 9, the condition of is also satisfied (Table 4), and the plot of versus in Figure 10(f) shows a peak at . These evidently show a multifractal geometry of the pore space. There are and for each model as well, indicating that with the increasing , the generalized fractal dimension decreases rapidly when and slowly for . Meanwhile, and are exhibited for the Hölder exponent, which is similar to the behaviour in the generalized fractal dimension. These features are followed by the general shape of the singularity spectrum, versus , with a narrow range for (the left branches of Figure 10(e) curves) compared to the range for (the right branches of Figure 10(e) curves), indicating a more significant heterogeneity of small pores than heterogeneity of large pores [58]. The increase in and with the worsen of the resolution demonstrates that the heterogeneity in large pores becomes stronger in upscaling (smoothing) of the voxels. However, the heterogeneity of small pores follows a more complex dependency on resolution, since values of and reach the lowest level when the resolution is 16 μm. The possible causes of this behaviour will be further investigated.

3.2.2. Multifractal Spectra of Flow Velocity

Results of the multifractal analysis on the simulated flow velocity were also obtained as generalized dimension spectra, , the Hölder exponent spectra, , and singularity spectra, . Typical data are listed in Table 5. Similar to that for the pore geometry, a turning point exists in the dependency of spectra on the resolution, as indicated by the lowest level of and also in M48. In considering of such a turning point, we divide the spectra into two groups and plot them, respectively, in Figures 11 and 12.

Curves of shown in Figures 11(a) and 12(a) for the flow velocity have similar shape in comparison with that for the pore geometry (Figure 10(a)), while exhibiting a significantly larger range and steeper change around . This is also indicated by the similar data of , , , and but significantly higher values of in Table 5, compared with that in Table 4. This difference also exists in the Hölder exponent spectra shown in Figures 11(c) and 12(c), comparing with Figure 10(c).

With regard to the group of M256 to M48, and of the flow velocity have negative correlations with the resolution, as indicated in Figures 11(b), while almost keeps constant. The values of , , and in this group are negatively correlated with resolution, as indicated in Figure 11(d). For the group of M64 to M16, all values of , , , and of the flow velocity are positively correlated with the resolution as exhibited in Figures 12(b) and 12(d). Both and in this group are almost constant.

The singularity spectra of the flow velocity shown in Figures 11(e) and 12(e) for different groups are similar, forming single-peak curves, which are generally like the curves in Figure 10(e) for the pore geometry. However, negative values of the Hausdorff dimension, , could be estimated when the Hölder exponent, , is high. As interpreted by Mandelbrot [59], the positive Hausdorff dimension defines a typical distribution of the measure, whereas the negative exists in the latent portion and refers to high sampling variability. The physical significance of the negative fractional dimension has been investigated well in studies of turbulence and diffusion-limited aggregation [59]. In this study, the pore geometry yields , while the flow velocity results in for some large values. Such a difference may indicate a stronger variability in the flow velocity. The curves of versus shown in Figures 11(f) and 12(f) for the flow velocity are also similar, forming a sharper peak at than that shown in Figure 10(f) for the pore geometry.

Data in Table 5 () and plots in Figures 11 and 12 indicate a multifractal distribution of the flow velocity simulated from all of the digital twin models. Both the generalized dimension and the Hölder exponent of the flow velocity vary in a larger range for than the range for , similar to that of the pore geometry. It indicates that the distribution of low-velocity nodes has stronger multifractal property than the distribution of high-velocity nodes. The nonlinear dependency of and values on the resolution has a turning point at M48 (16 μm) for the flow velocity, which is the same as the turning point for the pore geometry. It may imply that the change in the multifractal property of the pore geometry with the resolution worsening also leads to a change in the multifractal property of the flow velocity. However, when , the generalized dimension, , the Hölder exponent, , and the Hausdorff dimension, , of the flow velocity do not significantly change with the resolution, which is different from the behaviour of those for the pore geometry. This implies that the resolution-dependent multifractal geometry of large pores does not significantly affect the multifractal distribution of high-velocity flow.

4. Discussion

Resolution of a digital rock twin controls the identified porosity. In our study, both of the 3D porosity and 2D porosity fluctuate slightly when the resolution becomes coarse from 3 μm (M256) to 16 μm (M48) and decrease rapidly when the resolution changes from 16 μm (M48) to 48 μm (M16) (Figure 3). The difference in the porosity among models of M256, M192, M128, M96, M64, and M48 can be almost ignored, because it is less than 5% related to the average porosity. This small variation may be a result of the small change in the magnification () from M256 to M48. Among models of M48, M32, M24, and M16, the difference in porosity becomes much larger related to the average porosity, as the magnification becomes bigger (). The investigation by Jouini et al. [30] covers a ratio (=15) of magnification for carbonate samples, which also results in a more significant variation in the estimated porosity, but there are stochastic behaviours, and the correlation could be positive or negative. The sandstone sample in our study has a large number of small pores (Figure 1) that may not be recognized as pores during binarization at a coarse-resolution level (Figures 2(e)2(h)). Thus, the negative correlation between porosity and resolution is a reasonable result. A similar relationship is also reported for carbonate samples in Shah et al. [31] and sandstone samples in Guan et al. [32].

The resolution-dependent porosity may further trigger the change in the estimated permeability with resolution, because the porosity usually is a major control of the permeability. Figure 6 sufficiently indicates such an impact and highlights a strong amplification effect of the resolution dependence from the porosity to the permeability, where the 2D porosity mainly ranges between 0.2 and 0.5, while the slice permeability covers several orders of magnitude. This strong amplification effect has also been reported by Shah et al. [31]. However, the 3D porosity is not the dominant factor for the holistic permeability, as indicated in the opposite responses of them on the change in resolution (Figure 5). A comparable influence factor is the equivalent hydraulic radius of pores, as that estimated in Section 2.1 and shown in Figure 5(b). It has been well known according to Poiseuille’s law that the permeability is linearly proportional to both of the pore radius squared, , and the porosity [60, 61], as follows: where is the porosity. Thus, we obtained data of the 2D porosity and average value of pores on slices in each twin model and then estimated the value of for each slice. The relationship between the slice permeability and is exhibited with log-log plots in Figure 13(a), even though this equivalent permeability does not really represent that at the REV scale. A positive relationship is indicated by the data points. The linear correlation between and implies a power function that can be expressed as: where and . Note that , which is different from Equation (16) where . Data points in Figure 13(a) do not soundly fall in the vicinity of this correlation line, leading to a low determination coefficient ( is less than 0.3). This may be a consequence of the random change in pore-scale flow in a limited space that is smaller than REV. We also estimated the average values for a whole twin model and compared it with the holistic permeability, as shown in Figure 13(b). Apparently, the data points do not follow Equation (16) or Equation (17) but approximately exhibit a positive linear relationship. Thus, Poiseuille’s law should be carefully used for a multifractal porous material.

In general, the multifractal spectra of the pore geometry shown in Figure 10 are similar to previous studies on 2D images [30, 62] and 3D models [20, 63]. The curves of and for 2D images exhibited a smoother shape with conditions of and , in comparison to those obtained from 3D twin models, indicating a more homogeneous distribution of the 2D pore geometry. The singularity spectra, versus , of the 3D pore geometry is close to a hyperbolic curve (Figure 10(e)) with a peak point almost satisfying but not fully symmetrical on sides of the peak. The left and right tails of the curve in Figure 10(e) have different lengths, where the right tail is generally longer. This indicates a greater heterogeneity in small pores, because the long tail on the right reflects the greater range in when (dominated by small pores) as shown in Figure 10(c). This characteristic has also been observed in previous studies [20, 30]. For the impact of resolution on the multifractal spectra of the pore geometry, Figure 10(a) shows the lowering of the curve with the increasing of the resolution (coarsening voxels), especially for the part of . A similar response of the multifractal spectra on the resolution was also reported in Jouini et al. [30], where the curve of 2D images moves downward by using a smaller magnification (resolution changes from fine to coarse). Karimpouli and Tahmasebi [20] provided curves for the 3D pore geometry of rock samples identified from different resolution models, where increases for but almost keeps steady for when the resolution changes from coarse to fine. Our study reveals that may also be higher with finer resolution for . Thus, fine resolution twin models are required in capturing the actual pore geometry of rocks. For the sandstone sample in this study, the heterogeneity of small pores follows a more complex dependency on resolution, since values of and reach the lowest level when the resolution is 16 μm (Figure 10), which may be the exact turning point. The mechanism of such a turning point is still unclear but be possibly linked with the size distribution of solid particles in the sandstone.

The multifractal characteristics of fluid flow at the pore-scale are controlled by the multifractal pore geometry. In comparison with the multifractal pore geometry, the 3D distribution of experiences a stronger multifractal feature, as revealed by wider ranges in and (compare Figures 11 and 12 with Figure 10), and the ranges are more sensitive to the resolution. However, for the part of , the and curves of are lower and less sensitive to the resolution than the and curves of the pore geometry. It also implies that the preferential flow passageways (high velocity) have a fractal dimension that is smaller than the pore geometry, which may be linked with the significant increase in the flow velocity from the solid surface to the centre of a large pore channel (Figure 8). The singularity spectra of the flow velocity also result in approximately hyperbolic curves of as shown in Figures 11(e) and 12(e), while having more asymmetrical shape on sides of the peak than the singularity spectra of the pore geometry. Comparable singularity spectra have been obtained by Jiménez-Hornero et al. [28] for idealized porous media using 3D multifractal analysis on LBM simulation results of the pore-scale flow. Different features in our study are the longer tail on the right and wider range of covered by the singularity spectra. The resolution does not significantly influence the left part of the singularity spectra, whereas the right tail shows a complex response to the resolution change. As indicated, the heterogeneity in the estimated flow velocity is enhanced by using a much coarse-resolution twin model, especially for low values (with respect to the right tail of a singularity spectrum).

Results of the resolution-dependent multifractal characteristics in our study may be limited in the sandstone sample extracted from the Digital Rocks Portal, and further investigations are expected to cover a wide range of porous media. Synthetically generated digital samples with a self-similar structure can be also used to check the generic effect of the resolution on multifractal porous.

We should make a remark about the limitation of our investigation. The original digital twin model of the sandstone sample, M256, has a resolution of 3 μm. This image quality seems to be not good enough to obtain the relivable permeability for the sandstone sample, as evaluated from the criterion suggested by Saxena et al. [33]: the ratio of the dominate pore throat size to the voxel size should be at least 10 to achieve a steady value of permeability. According to the equivalent hydraulic radius of pores, , estimated for M256, the average value is ~13 μm, implying that the dominate pore throat size is quite smaller than 10 times of the resolution. Therefore, the resolution of M256 may be also too coarse, since a dozen of voxels were suggested to fill the pore throat [33, 40, 64]. However, as indicated in Figure 5(b), the estimated holistic permeability was not sensitive to the voxel size when the resolution is smaller than 8 μm. The permeability of this sandstone samples may not be dominated by small pore throats but highly depends on large pores that are more than 10 times of 3 μm. This is an open question for further investigations. Another limitation of this study is that the REV size of the sandstone sample is a bit larger than the model size. Further studies are expected to investigate the effect of REV on the multifractal characteristics of digital rock twins.

5. Conclusions

In this study, we built up nine 3D digital twin models with resolutions of 3 μm, 4 μm, 6 μm, 8 μm, 12 μm, 16 μm, 24 μm, 32 μm, and 48 μm for a 768 μm-length sandstone sample, from original images of 3 μm-resolution. Pore-scale water flow in these twin models is simulated using LBM, to obtain flow velocities in the 3D space of pores. Multifractal analyses of the pore geometry and the velocity magnitude distribution are both implemented with the sandbox method. According to the study results, conclusions can be summarized as follows: (1)All the 3D porosity, equivalent hydraulic radius of pores, and holistic permeability are not sensitive to the change in resolution when the resolution is fine enough. As the resolution continues to coarsen, the 3D porosity and equivalent pore radius decrease and increase, respectively. The holistic permeability is more sensitive to the change in the pore radius and then follows a positive response to the coarsening process. For each model, the equivalent slice permeability (slice is smaller than REV) shows a fluctuation pattern along the macroflow direction, which is similar to the fluctuation of 2D porosity but has significantly larger ranges and different extreme points(2)The distribution of the flow velocity is significantly heterogeneous with preferential flow passageways in the 3D space. When the resolution of the twin model changes from fine to coarse, the average velocity magnitude in pores and the Darcy velocity are almost steady first and then increase in a nonlinear manner. The increase in the Darcy velocity is much less than velocity magnitude in pores(3)Multifractal spectra indicate that the pore geometry is multifractal, and the multifractal parameters are dependent on the resolution of the digital twin. The generalized fractal dimension, , and the Hölder exponent, , decay rapidly with the increasing around . The ranges of and for are nonlinearly dependent on the resolution, with a lowest point at the resolution of 16 μm. The singularity spectra, versus , of the 3D pore geometry is close to a hyperbolic curve with a peak point that moves with the varying resolution. Fine resolution digital twins are required in capturing the actual pore geometry of rocks(4)The multifractal flow at the pore-scale is controlled by the multifractal pore geometry but may have different resolution-dependency behaviours. In comparison with the multifractal pore geometry, the flow velocity distribution has a stronger multifractal feature. The singularity spectra of the flow velocity also result in approximately hyperbolic curves of , while having more asymmetrical shape on sides of the peak than that of the pore geometry. The left part of the singularity spectra is insensitive to the resolution, whereas the right tail shows a complex response to the resolution change, indicating that the heterogeneity in the flow velocity may be overestimated by using a coarse-resolution twin model

Data Availability

All data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

This research was supported by the Key Research Program of the Institute of Geology & Geophysics, CAS (no. IGGCAS-201903) and the Second Tibetan Plateau Scientific Expedition and Research Program (STEP) (Grant No. 2019QZKK0904).