Representative Elementary Volume of Rock Using X-Ray Microcomputed Tomography: A New Statistical Approach

Rock heterogeneity is a key parameter influencing a range of rock properties such as fluid flow and geomechanical characteristics. The previously proposed statistical techniques were able to rank heterogeneity on a qualitative level to different extents; however, they need to select a threshold value for determination of representative elementary volumes (REV), which in turn makes the obtained REV subjective. In this study, an X-ray microcomputed tomography (μCT) technique was used to obtain images from different porous media. A new statistical technique was then used to compute REV, as a measure of heterogeneity, without the necessity of defining a threshold. The performance of the method was compared with other methods. It was shown that the calculated sum of the relative errors of the proposed method was lowest compared to the other statistical techniques for all tested porous media. The proposed method can be applied to different types of rocks for more accurate estimation of a REV.


Introduction
Accurate determination of rock heterogeneity is critical for a variety of industrial applications; for instance, it plays a key role in determining the reservoir's ability to recover oil and gas [1][2][3][4], carbon geostorage efficiency [5][6][7][8], contaminant mitigation and natural source zone depletion [9,10], water discharge and extraction rates [11,12], or geothermal energy production feasibility [13][14][15]. It is thus essential to understand rock heterogeneity in detail so that reliable predictions can be made and the targeted processes can be further optimised.
Heterogeneity is formed as a result of a verity of geological processes such as deposition, diagenesis, erosion, and structural deformation that ultimately control the geometry of sedimentary deposits [11]. Heterogeneity, thus, is a multiscale phenomenon that reflects the complexity of a geological formation. It is significant at different length scales: from submicron subgrain scale, which may include biogeochemical features such as skeleton structures of biological species forming limestone [16,17] or grain distribution effects [18,19], to micrometre heterogeneities induced by different grain shapes and sizes [20], to eventually millimetre to centimetre fractures in the rock generated by geomechanical processes [21,22] or at even larger scale faults, which have been created by tectonic movements [17,23].
Knowledge of flow in micrometre pore scale is required to understand macroscopic characteristics of flow in porous media; thus, prior to performing upscaling exercises of petrophysical properties [24][25][26][27], the accurate characterisation of the porous media at pore scale is crucial for quantifying the heterogeneity of a large system. The pore-scale properties can then be upscaled to core scale or even larger scales. Several works have been conducted to deal with upscaling nonadditive properties of heterogeneous porous media, such as permeability, from the pore scale to the core scale, for example, applying dual-scale pore network models linking pores at different length scales [28] or a direct numerical simulation method for coupling pores at different length scales [4,29].
One of the methods for quantifying the heterogeneity is using the representative elementary volume (REV), which has been applied to different rock types [30][31][32]. The REV is defined as the minimum volume of a rock that is representative of any larger volume, provided that the next level of heterogeneity at a larger scale has not been reached [33]. This concept is visualised in Figure 1 by plotting porosity against rock volume. The point at which the parameter, i.e., porosity in Figure 1, becomes constant identifies the REV. While a practical REV is required for performing reservoir simulations and providing geological and fluid dynamical predictions, due to experimental limitations originated partly from rock heterogeneity, a constant value cannot be generally obtained. Statistical methods such as direct observation analysis and nonlinear regression modellings, which employ standard deviation of porosity, coefficient of variation of porosity, or relative gradient error of porosity as predictors [20,[33][34][35], are therefore used to provide an estimate of REV. Nevertheless, these methods have some limitations in finding the value of REV. Visual observation, as discussed by Costanza-Robinson et al. [20], has the least accuracy, and the other three methods highly depend on the way of choosing their threshold values [34,35].
In this work, an X-ray microcomputed tomography (μCT) technique is employed to acquire high-resolution 3D images from different porous media for observing detailed pore morphology at micrometre to millimetre scale. These digital porous media were then subsampled, and their REVs were computed using different methods and compared with a new method proposed in this study. The proposed method uses linear regression to determine the volume at which the standard deviation of porosity is minimised. A weighted regression is then employed in which the estimated REV is imputed, and the corresponding porosity is estimated.

Methodology
2.1. Samples. Five porous media with different level of heterogeneities, as presented in Table 1, were selected including a highly homogeneous glass bead pack as a benchmark, a sand pack, two outcrop sandstones (Bentheimer and St. Bees), and one outcrop carbonate (Mount Gambier). In the case of the outcrop rocks, cylindrical core plugs of approximately 5 mm in diameter and 5-10 mm in length were cored from larger blocks. The unconsolidated materials were packed in hollow plastic cylinders.
2.2. X-Ray Microcomputed Tomography Imaging. The samples were imaged in a dry state with an Xradia Versa XRM500T μCT instrument at a resolution of~4 μm 3 . The subvolumes of the samples (~3 mm 3 ) were imaged to minimise beam hardening, i.e., to improve image quality [36]. The raw images were then filtered using a 3D nonlocal means filter [37] and segmented with a watershed algorithm [38], as shown in Figure 2. 2.3. Digital Measurements. The porosity was calculated by the ratio of the number of voxels in the void phase (pore space) to the total number of voxels (bulk volume) in the segmented image. As shown in Figure 1, the porosity (ϕ i ) can vary with sample volume (V) because of heterogeneity. The images were randomly subsampled into a number of different subvolumes (subsamples), with the condition of no overlap. Figures 3 and 4 exemplify the subsampling method and illustrate the pore space of the digital Bentheimer subsamples with different sizes. Table 2 provides the size of different subsamples. In total, 105 subsamples were acquired. The subsamples ranged from 100 3 voxels to 300 × 600 × 600 voxels, which corresponded to a volume range of~0.04 mm 3 tõ 10.38 mm 3 for the subsamples with different image resolutions (voxel sizes).

Results and Discussion
3.1. Visual Interpretation. One of the simplest techniques to determine REV is plotting porosity values against sample volume, as suggested by Bear [33]. Application of this method to the digital porous media is displayed in Figure 5. It is seen that (a) there is a notable variance of porosity associated with the samples from the smallest volume, which is related to heterogeneity [34], and (b), as expected, this variance narrows with the increase of volume. This implies that a REV can be estimated at sufficiently low variance. In terms of point (a), the variance is particularly significant for volumes smaller than 0.1 mm 3 , and it ranges from ±10% for the highly homogeneous glass bead pack to ±100% for the highly heterogeneous carbonate. The variance of porosity can be thus used as a qualitative indicator for rock heterogeneity.
The results of REV assessment, using Figure 5, by visual inspection, is provided in Table 3. The average porosity ( Figure 5) follows the trend described by Bear [33]. However, this is a highly subjective analysis. In particular, the technique of looking at the data points does not provide a clear indication of where exactly the REV need to be placed.

Standard Deviation Analysis.
Another method for identifying the cut-off (i.e., the REV value) is applying a mathematical criterion. The unbiased standard deviation s, Equation 2 Geofluids (1), [39] was used to perform the analysis on the samples.
Recalling that the sample volume was subsampled into the i th size, where s is the standard deviation (unbiased estimate), ϕ i is the porosity measured for the i th subvolume (i = 1, 2, 3, ⋯k), ϕ is the arithmetic average of the porosity of the i th subvolume, N is the total number of subsamples taken, i is the smallest subsample, and k is the largest subsample; the porosity variance as a function of volume can then be assessed by plotting s i versus subsample volume ( Figure 6).  Figure 2: From left to right: raw 2D image slice through each porous medium (dark = pore space; grey and white = solid phase), segmented 2D image slice (blue = pore space and black = solid phase), 3D pore spaces (blue = pore space and transparent = solid phase), and 3D rock matrix (different colours indicate different grain sizes).

Geofluids
As seen in Figure 6, the standard deviation rapidly drops with increasing subvolume, but the slope of the decay also rapidly decreases, approximately following a power-law relation s = AV b , where A and b are the least square nonlinear regression fitting parameters and V is the subsample volume (Table 4). While the power-law exponents b are all in the range (0:5 ± 0:2), coefficient A causes the largest difference in terms of REV. Note that the most homogenous sample, i.e., the glass bead pack, has the smallest coefficient and exponent, and the most heterogeneous sample, i.e., the carbonate, has the largest coefficient A. The exponent b is likely related to the grain size as the sand pack has the largest b and the largest grains.
The standard deviation curves in general give a good indication of heterogeneity, and different samples can be readily compared. However, a threshold value is required to determine the REV. Table 5 shows that a smaller threshold value resulted in a larger REV. At the same time, a smaller threshold had higher accuracy. REV thus strongly depends on the required accuracy. The experimental accuracy of a standard helium pycnometer, ±0.5% [40], was used as a threshold value, Figure 6 (dashed black line). It was found that the glass bead pack consistently had the smallest REV (0.4 mm 3 at 0.5% threshold) indicating that the sample was very homogenous, in contrast to Mt. Gambier, which always had the highest REV and thus was a very heterogeneous rock.
The results were consistent with data reported by Zhang et al. [35], who measured a REV of 5 mm 3 for Brent Triassic sandstone and a REV of 2 mm 3 for a crushed glass bead pack (0.5% threshold), and Stroeven et al. [39] conclude that the REV for a densely packed system is significantly lower than for a less dense pack, as presented in Figure 2 and Table 5. It was found that the grain size had likely a major impact on REV, which was also reflected in the relatively high REV values of the sand pack (Table 5); even though the sand pack was comparatively significantly more symmetrically formed, it had large grains.
3.3. Coefficient of Variation. The coefficient of variation, CV, Equation (2) [41], is unitless and therefore useful when comparing variation between datasets with substantial differences in their means. CV can be used for a REV analysis in a similar way as s.
The plot of CV as a function of subvolume ( Figure 7) shows that CV rapidly decreases with increasing sample volume. The CV values again can be approximated using powerlaw relations (Table 6), similar to s, as CV and s are closely related (Equations (1) and (2)). The exponents and the ratios between the coefficients of CV (Table 6) and s ( Table 4) fitting equations are similar, although the nominal values are different. Similar to the standard deviation analysis, a cutoff value needs to be selected to obtain REV. Three thresholds (Table 7) were tested, and it was again observed that REV strongly depended on the threshold value. Table 7 indicates that the carbonate rock had the largest heterogeneity as expected, followed by St. Bees sandstone, while the glass bead pack was the most homogeneous media. By looking at the μCT images, it seemed that the Bentheimer was more homogeneous than St. Bees because of its more symmetrical structure. However, the grain size in St. Bees was substantially smaller than in Bentheimer, as seen in Figure 2; thus, the grain size was a key factor, which strongly influenced REV.
Zhang et al. [35] measured a REV of 0.02 mm 3 for a Brent Triassic sandstone (porosity = 15:2%) at a CV = 0:2. Their result is consistent with the observations in this study; nevertheless, using CV = 0:05 delivers results that are more compatible with pycnometric measurements ( Table 7). The consequence of using such a CV cut-off is, however, noticeable; REV increases~100-fold, assuming the sandstones are somewhat similar, which seems likely. This threshold value can thus considerably change the REV and has to be selected carefully, and it needs to be compatible with any subsequent analysis, e.g., NMR measurements, which can observe individual atoms [42], and may need a different accuracy level than for instance capillary pressure-water saturation measurements where volume balances in core plugs are considered [43].

Relative Gradient Error
Criterion. An alternative method for determining REV is the relative gradient error (ε g ) analysis, Equation (3) [20]: where ϕ is the porosity, i is the subsample number, and V is the volume of the subsample. Figure 8 presents the plot of computed ε g for all samples versus sample volume; again, it is seen that ε g rapidly drops with increasing subvolume, and the data are much closer to   Geofluids each other especially after 1 mm 3 subvolume. This is also reflected in the statistical fits through the data points (Table 8), where generally high power-law exponents (~-0.5 to -2.2 with a median b value of -1.6) were calculated. Similar to the other statistical methods, to determine REV, a threshold value had to be selected, which similarly showed a significant influence of the threshold value (Table 9). This influence was, however, smaller than where s or CV analysis was used, mainly because of the lower sensitivity of ε g analysis, i.e., more similar numbers and curves generally need to be compared in the ε g analysis. For instance, the glass bead pack had ε g values smaller than 0.1 mm 3 for all thresholds tested while the REVs of the other four samples could hardly be accurately distinguished with the ε g approach. Note that Costanza-Robinson et al. [20] prescribed ε g < 0:2 as a threshold value, which, however, led to low accuracy in this case (Table 8). It, therefore, seems that the s and CV analysis are both superior to the ε g method while the relative gradient error method is still superior to the visual method, consistent with Costanza-Robinson et al. [20].

Regression Modelling.
To avoid the threshold choice in the above methods and therefore have a unique solution to REV value, a new method to estimate the volume of each rock is proposed in this study, at which point the variance of the porosity is minimised. To achieve this, a linear regression model is used, where volume is modelled as a function of the standard deviation of porosity. The standard deviations of the smallest subvolumes are identified as outliers and excluded from the model as they added unnecessary noise. To ensure the normality of the residual errors, the volume is modelled as log e ðVolumeÞ: where y is the volume, x is the standard deviation of the corresponding volume, β 0 is the constant, β 1 is the slope parameters, and ϵ is the error term of the regression equation.
To minimise the standard deviation in this equation, it is assumed to have a theoretical value of zero. When x is zero, the equation becomes log e y = β 0 , and therefore, the volume       7 Geofluids at which the standard deviation of the porosity is minimised is e β 0 , i.e., the exponential of the y-intercept. While this provides an estimated REV, the percent porosity of which the REV corresponds to is also calculated by employing a weighted linear regression of present porosity as a function of log e ðVolumeÞ: where p is the percent porosity, α 0 represent the constant, α 1 is the slope term of this model, and ϵ is the error term. The regression is weighted by the reciprocal of the standard deviation at each volume. Once the constant and slope terms are calculated, the REV for that particular rock is substituted into the equation, thus returning its corresponding percent porosity. The results for each rock are presented in Table 10 in their untransformed original units. The backtransformed regression equations are then shown visually in Figure 9.
3.6. Comparison of the Techniques. Sum of the relative errors (unitless ratio) of regression [44], obtained from each technique, is used to compare the above techniques appropriately: where abs is the absolute value, Re is the sum of the relative error of regression, f i is the model estimate, and Y i is the exact value. Table 11 lists the sum of relative errors obtained from each model fit. As seen in Table 11, the proposed regression technique presented the lowest errors of the regression. Note that the average porosity was used in the model to have only one measurement per subvolume to be consistent with other techniques. In addition, in the majority of cases, the proposed method predicted a REV larger than the other methods. However, the comparison of the estimates with the plotted data indicated that the proposed regression method was very sensible both in terms of prediction and error propagation. Many of the subsamples from the rocks still showed a considerable variation of the porosity at the estimates of REV obtained by other techniques. A notable advantage of the proposed method is that it removes the subjectivity of either visually estimating REV or choosing a threshold to estimate REV.

Conclusion
Quantification of rock heterogeneity can be accomplished using a REV, i.e., the larger the REV, the higher the heterogeneity. A range of porous media was imaged using the μCT technique and then subsampled (105 subsamples for each porous medium, a total of 525 subsamples). Different statistical methods with which REV can be estimated, including visual interpretation, standard deviation analysis, coefficient of variation analysis, and relative gradient error   9 Geofluids criterion analysis, were tested on the subsamples and compared with the new regression method proposed in this study.
The results indicate that the visual inspection approach is the least accurate approach. The relative gradient error, the standard deviation, and the coefficient of variation analyses have a higher degree of accuracy than the visual inspection approach. However, the results highly depend on the selected threshold values. Furthermore, it is not possible to determine a true REV in a theoretical sense or base on selecting threshold values. The proposed regres-sion modelling method, however, does not rely on a visual inspection and threshold value selection. The sum of relative errors of regression is also the lowest using the proposed technique. The method gives larger REV, which is satisfactory as many of the subsamples from different rocks show considerable variation in porosity at the estimated REV obtained by other techniques. In addition, it is shown that the grain size has a profound impact on REV, i.e., the larger the grains, the larger the REV, and the samples with a very ordered and symmetric structure can have a large REV if they contain large grains. Mt. Gambier carbonate ε g = 0:0100 V -1.573 0.0100 -1.573 Table 9: Influence of threshold value on REV; ε g analysis. Associated equations are shown in Table 8 and associated data in Figure 8.

Data Availability
The data supporting the conclusions of the study is already presented in the manuscript, and readers can download from the paper or can contact the authors.

Conflicts of Interest
There are no conflicts of interest.