Extended power-law scaling of heavy-tailed random air-permeability fields in fractured and sedimentary rocks

We analyze the scaling behaviors of two fieldscale log permeability data sets showing heavy-tailed frequency distributions in three and two spatial dimensions, respectively. One set consists of 1-m scale pneumatic packer test data from six vertical and inclined boreholes spanning a decameters scale block of unsaturated fractured tuffs near Superior, Arizona, the other of pneumatic minipermeameter data measured at a spacing of 15 cm along three horizontal transects on a 21 m long and 6 m high outcrop of the Upper Cretaceous Straight Cliffs Formation, including lowershoreface bioturbated and cross-bedded sandstone near Escalante, Utah. Order q sample structure functions of each data set scale as a power ξ(q) of separation scale or lag, s, over limited ranges of s. A procedure known as extended self-similarity (ESS) extends this range to all lags and yields a nonlinear (concave) functional relationship between ξ(q) andq. Whereas the literature tends to associate extended and nonlinear power-law scaling with multifractals or fractional Laplace motions, we have shown elsewhere that (a) ESS of data having a normal frequency distribution is theoretically consistent with (Gaussian) truncated (additive, self-affine, monofractal) fractional Brownian motion (tfBm), the latter being unique in predicting a breakdown in power-law scaling at small and large lags, and (b) nonlinear power-law scaling of data having either normal or heavy-tailed frequency distributions is consistent with samples from sub-Gaussian random fields or processes subordinated to tfBm or truncated fractional Gaussian noise (tfGn), stemming from lack of ergodicity which causes sample moments to scale differently than do their ensemble counterparts. Here we (i) demonstrate that the above two data sets are consistent with sub-Gaussian random fields subordinated to tfBm or tfGn and (ii) provide maximum likelihood estimates of parameters characterizing the corresponding Ĺ evy stable subordinators and tfBm or tfGn functions.


Introduction
Many earth and environmental (as well as physical, ecological, biological and financial) variables exhibit power-law scaling of the following type.Let (1) be an order q sample structure function of a random function Y (x) defined on a continuum of points x in one-or multidimensional space (or time), Y n (s) = Y (x n + s • m) − Y (x n ) being a sampled increment of Y (x) over a separation distance (lag) s in one or multiple directions, defined by one or more unit vectors m, between two points and N (s) the number of measured increments.Power-law scaling of Y (x) is described by S q N (s) ∝ s ξ (q)  (2) where the power or scaling exponent, ξ (q), is independent of s.When the scaling exponent is linearly proportional to q, ξ (q) = H q, Y (x) is interpreted to be a selfaffine (additive, monofractal) random field (or process) with Hurst exponent H .When ξ (q) varies nonlinearly with q, Y (x) has traditionally been taken to represent multiplicative, Published by Copernicus Publications on behalf of the European Geosciences Union.
Power-law scaling is typically assessed by employing the method of moments to analyze samples of measured variables.This entails inferring sample structure functions (Eq. 1) for a set q 1 , q 2 , ..., q n of q values at various lags.The structure function S q i N is related to s by linear regression on a log-log scale, the power ξ (q i ) (i = 1, 2, . . ., n) being set equal to the slope of the regression line.Linear or nearlinear dependence of log S q i N on log s is typically limited to intermediate ranges of separation scales, s I < s < s II , outside of which power-law scaling breaks down.The lower and upper limits, s I and s II respectively, which demarcate the range of power-law scaling are defined theoretically or, in most cases, empirically (Siena et al., 2012;Stumpf and Porter, 2012).Benzi et al. (1993a, b) provided empirical evidence that a procedure they had termed extended self-similarity (ESS) allows widening significantly the range of lags over which velocities in fully developed turbulence (where s I is taken to be governed by the Kolmogorov's dissipation scale) scale in a manner consistent with Eq. ( 2).Writing Eq. ( 2) as S n (s) = C(n) s ξ (n) and S m (s) = C(m) s ξ (m) , solving one of these equations for s and substituting into the other yields the ESS expression S n (s) ∝ S m (s) β(n,m)  (3) where β(n, m) = ξ(n)/ξ(m) is a ratio of scaling powers.Although the literature does not explain how and why Eq. ( 3) should apply to lags s < s I and s > s II where power-law scaling (Eq.2) breaks down, it nevertheless includes numerous examples demonstrating this to be the case.In addition to the classic case of turbulent velocities (Chakraborty et al., 2010), these examples include geographical (e.g.Earth and Mars topographic profiles), hydraulic (e.g.river morphology and sediment dynamics), atmospheric, astrophysical, (e.g.solar quiescent prominence, low-energy cosmic rays, cosmic microwave background radiation, turbulent boundary layers of the Earth's magnetosphere), biological (e.g.human heartbeat temporal dynamics), financial time series and ecological variables (see Guadagnini and Neuman (2011), Leonardis et al. (2012) and references therein).In virtually all these examples, ESS yields improved estimates of ξ (q) and shows it to vary in a nonlinear fashion with q, a finding commonly taken to imply that the variables are multifractal.Yet computational analyses by Guadagnini and Neuman (2011) have shown that this need not be the case: they found signals constructed from sub-Gaussian processes subordinated to truncated (additive, self-affine, monofractal) fractional Brownian motion (tfBm) to display ESS scaling as well as typical symptoms of multifractality, such as nonlinear scaling and intermittency, even though the signals differ from multifractals in a fundamental way (Neuman, 2010a(Neuman, , b, 2011;;Guadagnini et al., 2012).Siena et al. (2012) have pointed out that since multifractals and fractional Laplace motions do not capture observed breakdowns in power-law scaling at small and large lags, they cannot explain how and why ESS does.Instead, they have proven theoretically that ESS of data having a normal frequency distribution is theoretically consistent with tfBm.This allowed them to identify the functional form and estimate all parameters of the particular tfBm corresponding to log air-permeability data collected by Tidwell and Wilson (1999) on the faces of a laboratory-scale block of Topopah Spring tuff.In this paper we employ ESS to analyze the scaling behaviors of two log permeability data sets showing heavy-tailed frequency distributions in three and two spatial dimensions, respectively.One set consists of 1-m scale pneumatic packer test data from six vertical and inclined boreholes spanning a decameters-scale block of unsaturated fractured tuffs near Superior, Arizona (Guzman et al., 1996).Another set contains pneumatic minipermeameter data measured at a spacing of 15 cm along three horizontal transects on a 21 m long and 6 m high outcrop of the Upper Cretaceous Straight Cliffs Formation, including lower-shoreface bioturbated and cross-bedded sandstone near Escalante, Utah (Castle et al., 2004).Our analysis (a) demonstrates that the two data sets are statistically and theoretically consistent with sub-Gaussian random fields subordinated to tfBm or truncated fractional Gaussian noise (tfGn) and (b) provides maximum likelihood estimates of parameters characterizing the corresponding Lévy stable subordinators and tfBm or tfGn functions.

Theoretical background
We start by recounting the theory that underlies our analysis of the data.
It is possible to select a subordinator W 1/2 ≥ 0 having a heavy-tailed distribution other than Lévy such as, for example, a log-normal W 1/2 = e V with V = 0 and V 2 = (2 − α) 2 .Samples generated through subordination of truncated monofractal fBm in the above manner exhibit apparent multifractal scaling (Guadagnini et al., 2012).

Extended power-law scaling of sub-Gaussian processes subordinated to tfBm
It is important to note that whereas power-law scaling expressed by Eq. ( 2) implies ESS scaling in the form of Eq. ( 3), the reverse is not necessarily true because Eq. ( 3) follows from the more general relationship where f (s) is a (possibly nonlinear) function of s (Kozubowski and Molz, 2011;Siena et al., 2012).Following Neuman et al. (2012), we first consider subordinators W 1/2 ≥ 0 that have finite moments W q/2 of all orders q, such as the log-normal form mentioned earlier.Then, in a manner analogous to Siena et al. (2012), the central qthorder moments of absolute values of zero-mean stationary increments Y (x, s; λ l , λ u ) = W 1/2 G (x, s; λ l , λ u ) can be expressed as Here !! represents double factorial, i.e. q!! = q(q-2) (q-4). . . 2 if q is even and q!! = q(q-2) (q-4). . . 3 if q is odd, and The ratio between structure functions of order (q+1) and q is then where g (q) depends on the choice of subordinator but not on s.In the log-normal case where 8) and ( 9) that showing that log S q+1 is linear in log S q , in accord with the ESS expression Eq. ( 3), regardless of the choice of subordinator or the model employed for G (s; λ l , λ u ) 2 .On loglog plot, this line is characterized by a slope which tends to unity as q → ∞, being equal to 2 at q = 1.Equation ( 10) is a consequence of the equivalence between Eq. ( 8) and ESS expression Eq. ( 7) in which now f (s) = 2γ 2 (s; λ l , λ u ) .It shows that extended power-law scaling, or ESS, at all lags is an intrinsic property of sub-Gaussian processes subordinated to tfBm (or tfGn) with subordinators, such as the log normal, which have finite moments of all orders.
We noted earlier that, in the limits λ l → 0 and λ u → ∞, the TPV γ 2 i (s; λ l , λ u ) converges to a PV γ 2 i (s) = A i s 2H .It follows that Eq. ( 8) can be rewritten in terms of a power-law where it is clear that a log-log plot of S q versus s is linear at all lags and associated with a constant slope qH.
Following Neuman et al. (2012), we now consider subordinators W 1/2 ≥ 0 that have divergent ensemble moments W q/2 of all orders q ≥ 2α, as does the previously discussed Lévy subordinator with stability index α.In practical applications, Y (s; λ l , λ u ) q is typically estimated through a sample structure function where y m (x n , s; λ l , λ u ) denotes a collection of M < ∞ sets of N (s) < ∞ sampled increments each; for simplicity, we ignore possible variations of N (s) and x n with m.
Since order q ≥ 2α moments of w 1/2 m diverge while all moments of g m (x n , s; λ l , λ u ) converge, one can approximate Eq. ( 13) for a sufficiently large sample size N (s) by which, for finite M, is always finite.One can then write or, in analogy to Eq. ( 10), This indicates that S q+1 M (s; λ l , λ u ) on log-log scale, in accord with ESS expression (Eq.3), regardless of the functional form G (s; λ l , λ u ) 2 takes.The slope of this line is characterized by the same asymptotic behavior as that observed before.The approximate equivalence between Eq. ( 14) and the ESS expression (Eq.7), where f (s) = 2γ 2 i (s; λ l , λ u ) , is the basis for Eq. ( 16) and its asymptotic tendency.It follows that extended power-law scaling, or ESS, at all lags is an intrinsic property of samples from sub-Gaussian processes subordinated to tfBm (or tfGn) with subordinators, such as Lévy, which have divergent ensemble moments of orders q ≥ 2α.
Note that in the limits λ l → 0 and λ u → ∞, Eq. ( 14) becomes a power-law rendering log S q linear in log s with constant slope qH.

Analysis of log air permeabilities from borehole tests in unsaturated fractured tuff near Superior, Arizona
We analyze (natural) log air-permeability (Y = log k, k being permeability) data from unsaturated fractured tuff at a former University of Arizona research site near Superior, Arizona.Our analysis focuses on log k values obtained by Guzman et al. (1996) from steady state interpretations of 184 pneumatic injection tests in 1-m long intervals along 6 boreholes at the site (Fig. 1).Five of the boreholes (V2, W2a, X2, Y2, Z2) are 30 m long and one (Y3) has a length of 45 m; five (W2a, X2, Y2, Y3, Z2) are inclined at 45   Guzman et al. (1996).Riva et al. (2012) analyzed the probability distributions of 1-m scale (natural) log k measurements and their increments at the ALRS.They hypothesized that the data derive from a Lévy stable distribution, estimated the parameters of this distribution by three different methods and examined the degree to which each distribution estimate fits the data.Treating the data as a sample from a sub-Gaussian random field subordinated to tfBm via a Lévy stable subordinator and setting i = 1 in Eq. ( 5), they obtained the following ML estimates of the parameters characterizing the process: H = 0.33, σ 2 G = 4.05, λ l = 0.48 m and λ u = 9.98 m, thus yielding A = 0.67 by virtue of Eq. ( 9) in Riva et al. (2012).
Here we analyze structure functions and scaling of the same data using the ESS approach.We focus here on parameter estimates obtained by Riva et al. ( 2012) using a maximum likelihood (ML) approach applied to a log characteristic function of an α-stable variable, X; ϕ is a real-valued parameter; sign(ϕ) = 1, 0, −1 if ϕ > 0, = 0, < 0, respectively; α ∈ (0, 2] is stability or Lévy index; β∈[−1, 1] is skewness parameter; σ > 0 is scale or width parameter; and µ is shift or location parameter.The authors found Y = log k − log k to fit Eq. ( 18) with parameter estimates α = 2.0±0.00,σ = 1.42± 0.15 and μ = 0.00 ± 0.29.Note that it is difficult to estimate  β reliably when α ≈ 2 because, at α = 2, the distribution is insensitive to β. Figure 2a compares the frequency distribution of the data with their ML estimated probability density function and Fig. 2b depicts a corresponding Q-Q plot.The fits are ambiguous enough to suggest that their near-Gaussian appearance could in fact indicate a Lévy stable distribution with α just slightly smaller than 2. That this is likely the case follows from the tendency of α, fitted to the distributions of log k increments, to increase from 1.46 ± 0.21 at 1 m lag through 1.84 ± 0.16 at lag 2 m and 1.91 ± 0.12 at lag 3 m to 2 at lags equal to or exceeding 4 m.Increments corresponding to lags smaller than 4 m are thus clearly heavy tailed (and hence non-Gaussian) as evidenced further by Fig. 3, which compares frequency distributions and ML estimated probability density functions of Y = log k − log k data and log k increments at lags 1 m, 2 m and 5 m.Had the original log k data been genuinely Gaussian, the same would have to be true for their increments.Figure 4 depicts omnidirectional structure functions, S q N , of orders q = 1, 2, 3, 4, 5 computed for the same data according to Eq. ( 12).To compute them we ascribe each measurement to the midpoint of the corresponding 1-m scale borehole test interval.We then associate (as is common in geostatistical practice) data pairs separated by distances of 0.5-1.5 m with a lag of 1 m, those separated by distances of 1.5-2.5 m with a lag of 2 m, and so on up to the largest separation distances of 29.5-30.5 m, which we associate with a lag of 30 m. Figure 5 shows that the number of data pairs associated in this manner with each lag is largest at intermediate lags, causing log k increments to be comparatively undersampled at small and large lags.Such undersampling may explain in part why the structure functions in Fig. 4 scale differently with separation scale at small, intermediate and large lags.Standard moment analysis would entail fitting straight lines to these functions at intermediate lags by regression and considering their slopes to represent power-law exponents ξ (q) in Eq. ( 2).However, deciding what constitutes an appropriate range of intermediate lags for such analysis would, in the case of Fig. 4, be fraught with ambiguity.
We avoid this ambiguity by plotting in Fig. 6 S q N versus S q−1 N for 2 ≤ q ≤ 5 on log-log scale for the entire range of available lags.Also shown in Fig. 6 are linear regression fits to each of these relationships, the corresponding regression equations and coefficients of determination, R 2 .As the latter exceed 0.99 in all cases, we conclude with a high degree of confidence that S q N is a power β(q, q − 1) of S q−1 N for 2 ≤ q ≤ 5 at all lags, in accord with ESS expression (Eq.3).This power, given by the slopes of the regression lines in Fig. 6, decreases from 1.66 at q = 2 through 1.29 at q = 3 and 1.18 at q = 4 to 1.13 to q = 5, appearing to tend asymptotically toward 1 with increasing q.Considering S q N to vary as Hydrol.Earth Syst.Sci., 16, 3249-3260, 2012 www.hydrol-earth-syst-sci.net/16/3249/2012/ 2 3 ξ(q) 0 1 0 1 2 3 4 5 q Fig. 7. ξ (q) as a function of q (symbols) obtained via ESS based on ξ (1) = 0.56 computed for Arizona data by method of moments.Solid line has slope ξ (1) = 0.56 and dashed line slope H = 0.33 estimated for these data based on our theory, using maximum likelihood, by Riva et al. (2012).a power ξ (q) of s according to Eq. ( 2) at intermediate lags, as suggested by Fig. 4, allows expressing the power of S q N in Eq. (3) as β(q, q −1) = ξ(q)/ξ(q −1).Asymptotic tendency of β(q, q − 1) toward 1 then implies asymptotic tendency of ξ (q) toward a straight line.This commonly observed tendency, which the multifractal literature attributes to divergence of higher-order moments, is according to our theory (Neuman, 2010a;Guadagnini and Neuman, 2011) unrelated to such divergence, arising instead from the presence of an upper cutoff scale, λ u .
Figure 4 includes two vertical broken lines demarcating a mid-range of lags within which log S 1 N appears to be quite unambiguously linear in log s.Fitting a straight line to the corresponding data by regression yields ξ (1) = 0.56 with a high coefficient of determination, R 2 = 0.97.This, together with values of β(q, q − 1) = ξ(q)/ξ(q − 1) corresponding to 2 ≤ q ≤ 5 in Fig. 6, allows us to compute ξ (q) for this entire range of q values, as depicted in Fig. 7. Figure 7 also includes for reference one straight line having slope ξ (1) = 0.56 and another having slope H = 0.33, estimated for the same data by Riva et al. (2012).It is evident that ξ(q) in Fig. 7 is nonlinear concave in q in the range 2 ≤ q ≤ 5. Though such nonlinear scaling is typical of multifractals or fractional Laplace motions, we have demonstrated theoretically earlier that it is in fact consistent with a random field subordinated to tfBm via a heavy-tailed subordinator.

Analysis of nitrogen minipermeameter data from sandstone near Escalante, Utah
Castle et al. ( 2004) describe nitrogen minipermeameter measurements conducted on a flat, nearly vertical outcrop of Straight Cliffs Formation sandstones about 10 km northwest of Escalante, Utah.The outcrop, measuring approximately 21 m across and 6 m high, includes a lower bioturbated facies and an upper cross-bedded facies (Fig. 8).A total of 515 permeability measurements were taken in triplicate at a sample spacing of 15 cm along three horizontal transects (380 measurements) and four vertical profiles (135 measurements).Castle et al. (2004) found that whereas sample statistics of (natural) log permeability, log k, vary depending on which facies are considered, the frequency distributions of horizontal log k increments in the two facies are similar.Lu et al. (2002) used a fBm model to generate log k increments within a mix of distinct facies.They showed that, when data from different facies are analyzed jointly, the simulated log k increments exhibit an apparent non-Gaussian distribution.They concluded that observed Lévy-like behavior of sample probability distributions of permeability data can in some cases be an artifact of mixing data from disparate facies.Accordingly, Moltz et al. ( 2007) focused their analysis on increments along horizontal transects D and H (Fig. 8) within the lower bioturbated facies.They found the horizontal log k increments to be well represented by a fractional Laplace noise model.We note however that this same model would not have allowed them to characterize statistically the log k data themselves.
In this paper we analyze the frequency distributions and scaling of log k values and their horizontal increments (a) along transects D and H within the lower bioturbated facies and (b) jointly along transects D, H and X (Fig. 8) in the two facies.We also attempted to perform a similar analysis of log k values and their increments along the four vertical transects at the site but found the corresponding samples too small to yield meaningful statistics.
Transect H contains 133 data points, transect D 136 points and transect X 111 points.In a manner consistent with Riva et al. (2012), we analyze the frequency distribution of Y = log k − log k and use the computer code STABLE (Nolan, 1997(Nolan, , 2001) ) to obtain reliable ML estimates of stable densities.Figure 9a compares the frequency distribution of Y data from transects D and H on semi-logarithmic scale with a probability density function (pdf) fitted to it via ML.Treating the data as if they were Lévy stable yields ML parameter estimates α = 1.99 ± 0.05, σ = 0.28 ± 0.02, β = 0 and μ = 0.00 ± 0.05.As α ≈ 2, the distribution appears to be Gaussian.Yet Kolmogorov-Smirnov and Shapiro-Wilk tests reject the Gaussianity hypothesis at a 0.1 % significance level.The frequency distribution of Y data from all three horizontal transects D, H and X in Fig. 9b is positively skewed with ML parameter estimates α = 1.20 ± 0.12, β = 1, σ = 0.39 ± 0.04 and μ = 0.726 ± 0.07.We conclude that the two facies contain distinctly different log permeability populations Y .
Figure 10 compares frequency distributions and ML estimated probability density functions of log k increments along transects D and H, and jointly along transects D, H and X, at horizontal lags of 0.15 m, 0.45 m, 1.5 m and 4.5 m.Whereas at small lags the two distributions are similar (Fig. 10a, b), at larger lags the joint set from both facies exhibits heavier tails.Kolmogorov-Smirnov and Shapiro-Wilk tests generally reject the hypothesis that the increments, at any lag, are Gaussian at a 0.1 % significance level.A χ 2 test applied to horizontal increments along transects D and H at a lag of 0.15 m by Castle et al. (2004) has shown them to be Gaussian only at a 51 % confidence level.
As shown in Fig. 11, ML estimates α of the Lévy index of log permeability increments along transects D and H vary from 1.89 ± 0.13 at horizontal lag 0.15 m through 1.86 ± 0.14 at lag 0.3 m, 1.66 ± 0.18 at lag 0.45 m, 1.86 ± 0.14 at lag 0.6 m, 1.82 ± 0.16 at lag 0.75 m, 1.99 at lag 0.9 m to 2.00 at larger lags.Hence, the distributions of the increments have heavier tails at small rather than at larger lags.ML estimates α obtained from all three horizontal transects oscillate around 1.75 without any identifiable trend.ML estimates σ of the scale parameter in Fig. 11 increase monotonically with lag toward a constant asymptote of 0.32 for data along transects D and H and 0.44 for data along transects D, H and X.Both phenomena are consistent with the observation of Lu et al. (2002) that mixing data from the two facies may cause the tails of incremental frequency distributions to increase.
Results based on data sampled along transects D and H in the bioturbated sandstone facies are consistent with a sub-Gaussian random field subordinated to tfBm via a Lévy stable subordinator.The observed increase in α with lag is consistent with a version of such a field considered by Riva et al. (2012).Following their approach, Eq. ( 6) allows us to estimate the associated Hurst coefficient from the log-log slope of σ (s) in Fig. 11 at lags small enough to avoid the asymptote.This slope yields an estimate H = 0.13.From Eq. ( 6) it follows that, asymptotically, σ 2 G = 2 σ 2 where G (s; λ l , λ u ) is our tfBm.This, coupled with our ML estimates of σ for the log k − log k data, yields σ 2 G = 2×(0.28) 2 = 0.16.Having thus estimated H and σ 2 G , we are now in a position to estimate the remaining parameters of the TPV γ 2 G (s; λ l , λ u ) of G (s; λ l , λ u ) defined in Eq. ( 5).Setting i = 1 in Eq. ( 5), we obtain the following ML estimates of the cutoff scales, λ l ≈ 0.0 m and λ u = 16.97 m (with 95 % confidence limits 3.45 m and 30.47 m; setting i = 2 yields a less satisfactory fit, suggesting that i = 1 is a better choice), thus yielding A = 2.05 × 10 −2 by virtue of Eq. ( 9) in Riva et al. (2012).Our estimate of λ l is consistent with the small support scale of the minipermeameter.Our estimate of λ u is slightly smaller than the lengths of the D and H transects (on the order of 20 m), as expected from theory (Guadagnini et al., 2012).Figure 12 depicts experimental scale parameters and their theoretical equivalents based on the above ML estimates of σ 2 G , H , λ l and λ u .Dashed curves in the figure represent 95 % confidence limits of corresponding λ u estimates.
Results based on data sampled jointly along transects D, H and X in the bioturbated and cross-bedded sandstone facies are not fully consistent with our theory, which considers both Y and its increments to have symmetric distributions.As the distributions of the corresponding increments are in fact symmetric, it is possible to treat these increments as random field subordinated to truncated fractional Gaussian noise (tfGn) forming truncated sub-Gaussian fractional Lévy noise (tfLn) as discussed by Riva et al. (2012).Such processes are characterized by Lévy indices α that are independent of lag.Repeat- q th order sample structure functions q th order sample structure functions Fig. 13.Sample structure functions of order q = 1, 2, 3, 4, 5, 6 of Utah data from transects D and H (bioturbated sandstone).Light vertical broken lines demarcate mid-range of lags within which heavy inclined broken line, with slope taken to represent ξ (1), was fitted to S 1 N .
ing the above procedure, we obtain estimates H = 0.21, σ 2 G ≈ 0.34, λ l ≈ 0.0 m and λ u = 29.04m (with 95 % confidence limits 16.23 m and 41.85 m), thus yielding A = 4.32 × 10 −2 by virtue of Eq. ( 9) in Riva et al. (2012).Though this estimate of H exceeds that obtained previously on the basis of data from transects D and H alone, both are small and indicative of strong anti-persistence typical of log permeabilities in fractured and porous rocks worldwide (Neuman, 1990).
Figure 13 depicts sample structure functions of order q = 1, 2, 3, 4, 5, 6 for the data collected along transects D and H. Vertical lines demarcate the mid-range of lags within which a regression line, the slope of which was taken to represent ξ (1), had been fitted to S 1 N .The latter was found to be ξ (1) = 0.12 with coefficient of determination R 2 = 0.93.This value is only slightly smaller than that obtained earlier from the log-log slope of σ (s) in Fig. 11. Figure 14 shows loglog plots of S q N versus S q−1 N for 2 ≤ q ≤ 6 and corresponding linear regression fits.The fits are characterized by coefficients of determination, R 2 , two of which exceed 0.98 and three 0.99.The slope of the fitted lines decreases from 1.86 at q = 2 through 1.39 at q = 3, 1.25 at q = 4, and 1.18 at  for 2 ≤ q ≤ 6. Solid lines represent indicated regression fits.Linear regression equations and related regression coefficients (R 2 ) are also reported.q = 5 to 1.15 at q = 6, appearing to tend asymptotically toward 1 as expected.Adopting the above value of ξ (1) = 0.12 allows computing ξ (q) for 2 ≤ q ≤ 6 using the ESS relationship β(q, q − 1) = ξ(q)/ξ(q − 1).The results are plotted in Fig. 15 together with straight lines having slopes ξ (1) = 0.12 and H = 0.13.It is clear that ξ(q) is nonlinear concave in q within the range 2 ≤ q Though such nonlinear scaling is typical of multifractals or fractional Laplace motions, we have demonstrated theoretically earlier that it is in fact consistent with a random field subordinated to tfBm via a heavytailed subordinator.
Qualitatively similar results (details not given) are obtained from structure functions of order q computed jointly for horizontal increments along transects D, H and X in the two facies.Following the above procedure we obtain ξ (1) = 0.26, consistent with an analysis of σ (s) which yields H = 0.21.Applying ESS yields a nonlinear concave functional form for ξ(q) in Fig. 16, which also depicts for reference straight lines having slopes ξ (1) = 0.26 and H = 0.21.

Conclusions
Our analyses lead to the following major conclusions: 1. Extended power-law scaling, commonly known as extended self similarity or ESS, is an intrinsic property of sub-Gaussian random fields or processes subordinated to truncated fractional Brownian motion (tfBm) or truncated fractional Gaussian noise (tfGn).Such fields and processes are theoretically consistent with standard power-law scaling at intermediate lags and with ESS at all lags, including small and large lags at which powerlaw scaling breaks down.
2. Multifractals and fractional Laplace motions are theoretically consistent with standard power-law scaling at all lags.As such, they neither reproduce observed Hydrol.Earth Syst.Sci., 16, 3249-3260, 2012 www.hydrol-earth-syst-sci.net/16/3249/2012/ breakdown in power-law scaling at small and large lags nor explain how ESS extends power-law scaling to such lags.
3. 1-m scale pneumatic packer test data from unsaturated fractured tuffs near Superior, Arizona, and nitrogen minipermeameter data from bioturbated and crossbedded sandstones near Escalante, Utah and their increments show heavy-tailed frequency distributions that can be fitted with a high level of confidence to Lévy stable distributions.
4. Order q sample structure functions of each data set scale as a power ξ (q) of separation scale or lag, s, over limited ranges of s.ESS extends this range to all lags and yields a nonlinear concave functional relationship between ξ (q) and q.
5. The data sets we analyze are consistent with sub-Gaussian random fields subordinated to tfBm or to tfGn via Lévy stable subordinators.
6.This consistency allows estimating all tfBm or tfGn parameters (most notably the Hurst exponent and upper/lower cutoff scales) solely on the basis of the corresponding truncated power variograms.
7. The consistency further implies that nonlinear scaling of both data sets, manifested in a nonlinear concave relationship between their power-law exponents ξ (q) and q, is not an indication of multifractality but an artifact of sampling as explained theoretically by Neuman (2010a) and Guadagnini et al. (2012).

Fig. 2 .
Fig. 2. (a) Frequency distribution (symbols) and ML estimated probability density function (solid curve) of Arizona data; (b) Q-Q plot of empirical data versus theoretical estimate of stable distribution.

Fig. 3 .
Fig. 3. Frequency distributions (symbols) and ML estimated probability density functions (curves) of Arizona Y = log k − log k data (red) and log k increments at lags s = 1 m (black), 2 m (green), and 5 m (blue).

Fig. 14 .
Fig. 14.Log-log variations of S q N of Utah data from transects D and H (bioturbated sandstone) with S q−1 N Figure15