Interactive comment on “ Extended power-law scaling of air permeabilities measured on a block of tuff ”

This paper represents a substantial contribution to the understanding of the oftenobserved power-law scaling of porous media (soil, aquifers, . . .) properties and specifically hydrogeologic variables such as permeability and porosity. The authors used the methods of moments and of the so called extended power-law scaling to identify power law scaling of log air permeability data collected (using four measurement support scales) on the faces of a laboratory-scale cube of tuff and published by Tidwell and Wilson in 1999. The methodology used for analyzing the data are presented in a clearcut manner and the discussion/conclusions are original bringing a new insights into the derivation and significance of these power-law scaling (and range of occurrence)


Introduction
The literature indicates (Neuman and Di Federico, 2003) that hydrogeologic variables exhibit isotropic and directional dependencies on scales of measurement (data support), observation (extent of phenomena such as a dispersing plume), sampling window (domain of investigation), spatial correlation (structural coherence), and spatial resolution (descriptive detail).Attempts to explain such scale dependencies have focused in part on observed and/or hypothesized powerlaw behaviors of structure functions of variables such as hydraulic (or log hydraulic) conductivity (e.g.Painter, 1996;Liu and Molz, 1997a,b;Tennekoon et al., 2003), spacetime infiltration (Meng et al., 2006), soil properties (Caniego et al., 2005;Zeleke andSi, 2006, 2007), electrical resistance, natural gamma ray and spontaneous potential (Yang et al., 2009) and sediment transport data (Ganti et al., 2009).Power-law behavior means that a sample structure function (1) Published by Copernicus Publications on behalf of the European Geosciences Union.
M. Siena et al.: Extended power-law scaling of air permeabilities measured on a block of tuff of order-q (for simplicity we limit our mathematical exposition to one dimension and our analysis of data to nonnegative values of q) scales according to S q N (s) ∝ s ξ (q)  (2) where Y (x) is the variable of interest (assumed to be defined on a continuum of points x in space or time), Y n (s) is a measured increment Y (s) = Y (x + s) − Y (x) of the variable over a separation distance (lag) s between two points on the x-axis, and N (s) is the number of measured increments.
When the scaling exponent (power) ξ(q) varies linearly with q, Y (x) is interpreted to form a self-affine (mono-fractal) random field and the slope H of the corresponding line is termed Hurst exponent.When the scaling exponent ξ(q) is a nonlinear function of q, Y (x) has traditionally been taken to form a multifractal field.A semi-empirical "universal" multifractal model due to Schertzer and Lovejoy (1987) relates ξ(q) to the Hurst exponent via H = ξ(1), as explained and illustrated by Seuront et al. (1999); some approximate H by dξ/dq near q = 0. Neuman (2010aNeuman ( , 2011) ) has shown theoretically and Neuman (2010b) and Guadagnini et al. (2011) have demonstrated numerically that signals derived from additive processes subordinated to a truncated version (tfBm) of additive, selfaffine fractional Brownian motion (fBm) scale in a manner similar to multifractals even as they differ from such multiplicative constructs in a fundamental way.Truncation is caused by lower and upper cutoff scales, proportional to data support (measurement scale) or resolution and domain size (sampling window scale), respectively.Their work suggests that nonlinear variations in ξ(q) with q need not represent multifractal scaling but could instead be an artifact of sampling from tfBm or fields subordinated to tfBm.
Power-law scaling is typically inferred from measured values of earth and environmental variables by the method of moments (M).This consists of calculating sample structure functions Eq. (1) for a finite sequence, q 1 , q 2 , ..., q n , of q values and for various separation lags.For each order q i the logarithm of S q i N is related to log s by linear regression and the power ξ(q i ) set equal to the slope of the regression line.Linear or near-linear variation of log S q i N with log s is typically limited to intermediate ranges of separation scales, s I < s < s II , where s I and s II are theoretical or empirical lower and upper limits, respectively.Breakdown in power-law scaling is attributed in the literature to noise at lags smaller than s I and to undersampling at lags larger than s II (Tessier et al., 1993).Yet noise-free signals subordinated to tfBm generated by Neuman (2010b) and Guadagnini et al. (2011) show power-law breakdown at small and large lags even when sample sizes are large.This breakdown is caused by cutoffs which truncate the fields at small lags proportional to the measurement and/or resolution scale of the data, and at large lags proportional to the size of the sampling domain, regardless of noise or undersampling.Though nonlinear variation of ξ(q) with q is also reproduced by the fractional Laplace model of Meerschaert et al. (2004) (see Kozubowski et al., 2006;Ganti et al., 2009), the latter does not include cutoffs and thus fails to reproduce observed breakdown in power-law scaling at small and large lags.Benzi et al. (1993aBenzi et al. ( ,b, 1996) ) discovered empirically that the range s I < s < s II of separation scales over which velocities in fully developed turbulence (where Kolmogorov's dissipation scale is assumed to control s I ) scale according to Eq. ( 2) can be enlarged significantly, at both small and large lags, through a procedure they called Extended Self-Similarity (ESS).ESS arises from the observation that structure functions of different orders, n and m, computed for the same separation lag are related by where β(n, m) = ξ(n)/ξ(m) is a ratio of scaling exponents.Benzi et al. (1996) introduced, and Nikora and Goring ( 2001) employed, a generalized form of ESS (G-ESS) according to which where The exponent ρ(p, q, n) is a ratio between deviations of structure functions of order p and q, respectively, from linear (monofractal or self-affine) scaling.Chakraborty et al. (2010) cite the success of ESS in extending observed scaling ranges, and thus allowing more accurate empirical determinations of the functional exponent ξ(q) for turbulent velocities.ESS has been reported to achieve similar results for diffusion-limited aggregates, natural images, kinetic surface roughening, fluvial turbulence, sand wave dynamics, Martian topography, river morphometry, gravel-bed mobility and atmospheric barometric pressure, low-energy cosmic rays, cosmic microwave background radiation, metalinsulator transition, irregularities in human heartbeat time series, turbulence in edge magnetized plasma of fusion devices and turbulent boundary layers of the Earth's magnetosphere (Guadagnini and Neuman, 2011).In all cases, ESS has revealed nonlinear variation of ξ(q) with q.Whereas the literature has interpreted this to imply that ESS applies to multifractals, Guadagnini and Neuman have shown that Eq. (3) works equally well when applied to signals derived from additive processes subordinated to tfBm.As the latter are not multifractal, neither must be processes revealed by ESS (or any other method of analysis) to yield nonlinear variations in ξ(q) with q.
In this paper we use three methods to identify power-law scaling of log air permeability data collected by Tidwell and Wilson (1999) on the faces of a laboratory-scale cube of Topopah Spring tuff: method of moments (M) and extended power-law scaling via ESS and G-ESS.Most published analyses of extended power-law scaling concern time series or one-dimensional transects of spatial data associated with a unique measurement (support) scale.We use instead data measured on diverse support scales and distributed in two or three dimensions across several faces of the cube.Our aim is to infer the scaling behavior of these data using all three methods, compare results among the methods and explore the dependence of corresponding scaling exponents on support scales and direction.
"In spite of several attempts to explain the success of ESS" cited by Chakraborty et al. (2010) the authors note that "the latter is still not fully understood and we do not know how much we can trust scaling exponents derived by ESS.It would be nice to have at least one instance for which ESS not only works, but does so for reasons we can rationally understand."Chakraborty et al. (2010) provide such a theoretical reason in the special context of one-dimensional Burgers equation.In contrast, they consider "the multifractal description of turbulence," with which ESS is commonly associated, to be "quite heuristic and arbitrary."Kozubowski and Molz (2011) note that Eq. ( 3) is obtained from Eq. ( 2) simply upon rewriting the latter as S n (s) = C(n)s ξ(n) and S m (s) = C(m)s ξ(m) , solving the first of these expressions for s and substituting into the second.Kozubowski and Molz point out further that whereas Eq. (2) implies Eq. (3) the reverse is generally not true, Eq. (3) being equivalent instead to S q (s) ∝ f (s) ξ (q)  (6) where f (s) is some, possibly nonlinear, function of s.This is seen upon rewriting Eq. ( 6) as S n (s) = C(n)f (s) ξ(n) and S m (s) = C(m)f (s) ξ(m) , solving the first for f (s) and substituting into the second.
After showing that our data behave as a sample from tfBm (a truncated self-affine process) we demonstrate in Appendix A that this process is consistent with Eq. ( 6) at all separation scales (lags s) and with Eq. (2) at intermediate scales (s I < s < s II ), as are most of our data.We thus explain why and how ESS works for our data at all scales.At intermediate scales where our data are consistent with both Eqs. ( 6) and (2), the definition of β(n, m) in Eq. (3) applies and allows us to compute ξ(n) uniquely and unambiguously upon obtaining ξ(1) independently from the data by the method of moments.Similarly, the definition of ρ(p, q, n) in Eq. ( 5) applies and allows us to compute ξ(n) uniquely and unambiguously for any n upon obtaining two of its values, ξ(p) and ξ(q) for some p and q, independently from the data by the method of moments.On the other hand, since our data are consistent with Eq. ( 6) but not with Eq. (2) at small (s < s I ) and large (s > s II ) scales, they are inconsistent with multifractals or fractional Laplace motions (Meerschaert et al., 2004;Kozubowski et al., 2006;Ganti et al., 2009) which theoretically scale according to Eq. ( 2) at all lags.In other words our data, being consistent with a truncated self-affine process, exhibit apparent rather than actual multifractal scaling at intermediate lags.The same likely holds true for other Gaussian or heavy-tailed earth and environmental variables (such as those listed earlier) that scale according to Eq. ( 2) at intermediate lags and according to Eq. ( 3) over an extended range of lags, a possibility noted earlier by Guadagnini and Neuman (2011).
2 Previous analyses of experimental data Tidwell and Wilson (1999) measured air permeabilities, k, on six faces of a block of Topopah Spring tuff (Fig. 1), extending 30 cm on each side, with the aid of a Multisupport Permeameter (MSP).Measurements were conducted at intervals of = 0.85 cm on a grid of 36 × 36 points along each face using four tip-seal sizes having inner radii r i = 0.15, 0.31, 0.63, 1.27 cm and outer radii 2 r i .As the precise nature and size of the support volume associated with each measurement is the subject of debate (Goggin et al., 1988;Molz et al., 2003;Tartakovsky et al., 2000;Neuman and Di Federico, 2003), we consider the inner radius of the tip-seal to represent a nominal measurement scale (data support) as proposed by Tidwell and Wilson (1999).We conclude from their analysis that measurements on face 6 of the block are less reliable than the rest and therefore limit our analysis below to those on faces 1-5.
Measured (natural) log permeability values, Y = ln k, were found to have bi-modal frequency distributions particularly at larger tip sizes (Fig. 2 of Tidwell and Wilson, 1999).This was deemed by them to be consistent with the geologic structure of the tuff sample within which regions of high (associated with pumice fragments) and low (corresponding to solid matrix) permeability could be visually identified.Tidwell and Wilson were able to fit spherical models with nuggets to sample variograms on all faces of the cube for each tip radius.The variograms were found to be isotropic in the xy plane of Cartesian coordinates on face 1 of the cube but anisotropic in the xz and yz planes on faces 2-5, with estimated ranges in the z direction about one half of those in the x and y directions.Sill and range estimates decreased and increased, respectively, with tip seal inner radius.For additional details the reader is referred to the above authors.

Identification of power-law scaling
To evaluate sample structure functions for the experimental data of Tidwell and Wilson (1999) according to Eq. ( 1) we compute directional increments, Y , of Y = ln k at various separation lags (taken to be integer multiples of grid spacing, , for each tip size) parallel to the x, y and z-coordinates on the faces of the cube.Figure 2 depicts variations in Y associated with lag s x = 8.5 cm along selected transects in the x direction on face 1 associated with the smallest and largest tip radii, r i = 0.15 and 1.27 cm.Clearly, increasing the tip radius results in smoother and more persistent variability of the increments.Figure 3 shows frequency distributions of similar increments along all x-directional transects on multiple faces and Maximum Likelihood (ML) fits of Gaussian probability density functions (pdfs) to these distributions.Both frequency distributions are symmetric about zero with humps reflecting the bi-modal distributions of Y identified by Tidwell and Wilson.Increments corresponding to the smallest tip size appear to be closer to Gaussian than those corresponding to the largest tip size, consistent with their finding that Y becomes increasingly bimodal with tip size.
r i = 0.15 cm : sample ML r i = 1.27 cm : sample ML

Analysis of face 1 data by method of moments
As lag increases the number N (s) of incremental data along all transects of face 1 decreases from 1260 corresponding to s x = s y = = 0.85 cm to 36 corresponding to s x = s y = 35 × = 29.75cm. Figure 4 depicts log S q N (s x ) as functions of log s x along all transects for 0.1 ≤ q ≤ 2.5 at each tip size.To identify a middle range of lags within which these relationships are linear we have fitted regression lines to the data within several such ranges and adopted those that yield the highest coefficients of determination for each tip size.These ranges, identified in Fig. 4 by dashed vertical lines, are on the order of (s I = 2 ) ≤ s x ≤ (s II = 6 ) or 1.7 cm ≤ s x ≤ 5.1 cm.The corresponding nonlinear variation of ξ(q) with q for the largest tip size (based on data such as those in Fig. 4d) is depicted in Fig. 5.The solid line has slope dξ/dq| q=0 = 0.74 which, if taken to represent a Hurst exponent H , implies a persistent signal consistent with that depicted in Fig. 2b.Values of ξ(q) in Fig. 5 start deviating from this solid line at about q ≈ 0.6 to become asymptotically linear in q at about q ≥ 3.5, as evidenced by the dotted line obtained through regression against these values.Results for other tip sizes and in the y direction (not reported) are qualitatively similar.Though such behavior would typically be interpreted to imply that increments of ln k are multifractal, we note that qualitatively similar scaling has been produced synthetically by Guadagnini and Neuman (2011, their Fig. 4) with a model in which Y is subordinated to tfBm, a truncated version of self-affine (monofractal) fBm.

Analysis of face 1 data by extended power-law scaling
Replotting the data in Fig. 4d (corresponding to the largest tip size) as log S q N versus log S q−1 N for 2.0 ≤ q ≤ 5.0 (at intervals of 0.5) reveals much less ambiguous power-law scaling over a much wider range of lags in Fig. 6.Equations of corresponding curves (regression lines on log-log scale) included in the figure are characterized by coefficients of determination, R 2 , that exceed 0.98 at all lags.Results of similar quality (not reported) have been obtained for all tip sizes and directions.The slopes of the regression curves, representing β(q, q −1) in Eq. ( 3), decrease asymptotically with q toward unity consistently with the asymptotic tendency of ξ(q) in ξ (q) y = 0.74 x y = 0.09 x + 0.78 Fig. 5. ξ(q) versus q evaluated for r i = 1.27 cm on face 1 along xaxis.Continuous line has slope similar to ξ(q) near q = 0. Dashed line has slope similar to ξ(q) for q ≥ 3.5.
Fig. 5 toward linear variation with q.In Appendix A we explain this behavior theoretically by demonstrating that tfBm scales according to Eq. (2) at intermediate lags and according to Eq. ( 3) at all lags.The fact that our data scale according to Eq. ( 2 1.E+2 q = 2.5 q = 3.0 q = 3.5 q = 4.0 q = 4.5 Fig. 6. S q N versus S q−1 N for 2.0 ≤ q ≤ 5.0 and r i = 1.27 cm evaluated on face 1 along x-axis.Linear regression equations and relative regression coefficients (R 2 ) are also reported.
As ξ(q) in Fig. 5 starts deviating from the solid line at approximately q ≈0.6 we set n = 0.5 in Eq. ( 5) and plot in Fig. 7 log G n,q+1 (s x ) versus log G n,q (s x ) for q = 1.0, 1.5, 2.0, ..., 4.0 corresponding to the increments in Fig. 4d.Included in Fig. 7 are equations of curves fitted to these log-log relationships by linear regression and associated R 2 values.The figure reveals extended power-law scaling with R 2 ≥ 0.98 over virtually the entire range of lags.As in the earlier case of β(q, q − 1), the scaling ratios ρ(q + q, q, n) diminish asymptotically toward unity as q increases.Similar behavior is observed in the case of other tip sizes.Resulting values of ξ(q) corresponding to the x and y directions on face 1, computed in a manner analogous to that described in the previous paragraph and identified as G-ESS, are plotted versus q in Fig. 8.
Figure 8 juxtaposes values of ξ(q) as functions of q within the range 0 ≤ q ≤ 5.0, in the x and y directions of face 1, obtained for the largest tip size by the method of moments and two methods of extended power-law scaling.We saw earlier that the latter two methods are much less ambiguous than the first in helping one to identify and quantify powerlaw scaling of structure functions at various orders q.As ESS requires only one reference value, ξ (1), to compute ξ(q) on the basis of β(q, q − 1) for any order q while G-ESS requires two such reference values, we consider the former more reliable than the latter.
ξ Fig. 8. ξ(q) versus q evaluated for r i = 1.27 cm on face 1 in x and y directions.
Figure 9 shows that values of β(q, q − q) and ρ(q + q, q, n) are relatively insensitive to tip size and direction.The same is not true for the scaling exponent ξ(q) which, as shown in Fig. 10, increases consistently with tip size.Though these results correspond to the x direction on face 1, they do not differ qualitatively from those corresponding to x and y on all five faces.This behavior translates into a consistent increase in the Hurst exponent H with tip size (from H = 0.13 for r i = 0.15 cm to H = 0.74 for r i = 1.27 cm), implying that averaging over larger and larger support volumes smoothes a signal and renders it more persistent.
r i = 0.15 cm r i = 0.31 cm r i = 0.63 cm r i = 1.27 cm ξ (q) Fig. 10.ESS estimates of ξ(q) versus q evaluated on face 1 in x direction and various r i .
for each tip size, yielding 12 sets of three directional increments for 4 tip sizes.Figure 11 depicts log-log plots of S q N (s z ) versus separation distance, s z , along the z direction for 0.1 ≤ q ≤ 2.5 corresponding to each tip size.In neither plot is it possible to identify an intermediate range of powerlaw scaling, most likely due to the reduced range of the increments in this direction.We therefore omit incremental data in the z direction from further consideration in this paper.Structure functions in the x and y directions (not shown) display behaviors qualitatively similar to those noted earlier in the x direction on face 1 (Fig. 4). Figure 12 compares values of ξ(q) obtained by each method on all available data with xdirectional values obtained on face 1 via ESS.Whereas face 1 values in Fig. 8 show no significant difference between directions x and y, the multiface values in Fig. 12 do suggest a slight directional dependence revealed, most likely, by the relatively large size of this sample.In general, multiface values of ξ(q) in Fig. 12 lie below the face 1 values in Fig. 8, reflecting the impact of sample size on the quantification of power-law scaling.
As in the case of face 1, multiface values of β(q, q − q) and ρ(q + q, q, n) are seen in Fig. 13 to be relatively insensitive to tip size and direction.
We note that there is no conflict between the ability of Tidwell and Wilson (1999) to characterize Y for any tip size by means of a stationary variogram and our finding that order-q structure functions of Y exhibit power-law scaling at intermediate lags with exponents that vary in a nonlinear fashion with q.Instead both, coupled with our finding that increments of Y associated with small tip sizes are approximately Gaussian, are consistent with a view of Y as a sample from tfBm (implying that Y is not multifractal).Such a sample is characterized by a truncated power variogram (Di Federico and Neuman, 1997) which is difficult to distinguish from stationary variogram models (Neuman et al., 2008) and exhibits power-law scaling with exponents that are nonlinear in q at intermediate lags (Neuman, 2010a(Neuman, ,b, 2011;;Guadagnini et al., 2011).The former implies that Gaussian samples commonly characterized in the literature by stationary variograms may in fact represent truncated self-affine fields, the latter implies that such samples may in turn display apparent multifractality as does Y in this paper.

Model identification and parameter estimation
As shown by Eq. (A2) in Appendix A the q-th order structure function of tfBm is completely defined by q and the ensemble (theoretical) truncated power variogram (TPV) γ 2 (s; λ l , λ u ).Since increments of Y associated with the smallest tip size have a near-Gaussian distribution, one should be able to estimate the parameters of this variogram by fitting such theoretical structure functions to their sample counterparts, S q N , for r i = 0.15 cm.As the number of data N needed to obtain stable S q N values increases with q, we limit our estimation of parameters to structure functions S 1 N and S 2 N of orders q = 1, 2.
Equations ( A4)-(A7) in Appendix A make clear that a TPV is defined by four parameters: the Hurst exponent H , coefficient A, upper cutoff λ u and lower cutoff λ l .We found Figure 11 r i = 0.15 cm r i = 0.31 cm Fig. 12. ξ(q) versus q evaluated for r i = 1.27 cm using all the available data in x and y directions.ESS estimates obtained in x direction on face 1 are included for comparison.
earlier from the slope of ξ(q) at small q that, for incremental data on face 1, H = 0.13 in the x direction and H = 0.09 in the y direction while, for incremental data on multiple faces, H = 0.08 in the x direction and H = 0.09 in the y direction.
All four parameters are linked by the relationship where σ 2 is the sill (asymptotic plateau) of the variogram γ 2 (s; λ l , λ u ).We expect the estimate of this sill to not differ significantly from the sample variance of the Y data.
We estimate parameters on the basis of sample variograms, γ 2 * (s; λ l , λ u ), of Y data computed from sample structure functions of first and/or second order, respectively, via according to Eq. (A2).We estimate the sill, σ 2 , by averaging values of γ 2 * (s; λ l , λ u ) corresponding to large lags, s, obtained from S 1 N and/or S 2 N in this manner.We then estimate the cutoffs (and A) through a maximum likelihood (ML) fit of γ 2 (s; λ l , λ u ) to γ 2 * (s; λ l , λ u ) where the first is a TPV model based either on Gaussian or on exponential modes as defined in Eqs.(A4)-(A7).The ML procedure consists of minimizing the log likelihood criterion (Carrera and Neuman, 1986) with respect to λ u and λ l subject to Eq. ( 7).Here γ 2 and γ * 2 are vectors of n discrete γ 2 (s; λ l , λ u ) and γ 2 * (s; λ l , λ u ) (b) (a) ρ (q+∆q,q,n) β (q,q-∆q) Figure 13 x-axis: r i = 0.15 cm r i = 0.31 cm r i = 0.63 cm r i = 1.27 cm y-axis: r i = 0.15 cm r i = 0.31 cm Fig. 13.(a) β(q, q − q) and (b) ρ(q + q, q, n = 0.5) versus q evaluated using all available data in x and y directions and various r i .
values, respectively, T denotes transpose, C γ = σ 2 γ V where C γ is the covariance matrix of errors in γ * 2 (resulting from log permeability measurement errors), σ 2 γ is estimated during inversion according to where J min is the minimum of J , and Vis a known symmetric positive-definite matrix.For simplicity we take errors in γ * 2 to be uncorrelated and set V equal to the identity matrix.
Applying the above procedure to face 1 yields a sill of 4.48 based on S 1 N as well as on S 2 N of increments parallel to the xaxis, 3.64 based on S 1 N and 3.88 on S 2 N of increments parallel to the y-axis, the variance of the corresponding Y data being 4.06.Applying the above procedure jointly to faces 1, 2 and 4 (where incremental data in the x direction are available, see Fig. 1) yields a sill of 3.52 based on S 1 N and 3.77 based on S 2 N of increments parallel to the x-axis, the corresponding variance of Y being 3.77; applying it to faces 1, 3 and 5 (where incremental data in the y direction are available, see Fig. 1) yields 3.52 based on S 1 N and 3.79 on S 2 N of increments parallel to the y-axis, the corresponding variance of Y being 3.91 while that of all Y data on faces 1-5 is 3.79.We conclude that to obtain consistent estimates of σ 2 it is best to consider jointly all data from faces 1-5 as we do below.
Due to the irregular behavior of S 1 N (s) and S 2 N (s) at large lags (Fig. 4a) we limit our ML estimation of cutoffs (and A) to lags in the range ≤ s ≤ 13 so that n = 13.Table 1 lists parameter estimates and corresponding 95 % confidence intervals for TPV models consisting of Gaussian and exponential modes obtained on the basis of S 1 N , S 2 N and both with xand y-directional increments.The table also lists J min , NLL, the determinant |Q| of the covariance matrix Q of λ u and λ l estimation errors, and the Bayesian model discrimination criterion KIC (Kashyap, 1982) where M = 2 is the number of parameters.Values of quantities obtained on the basis of S 1 N , S 2 N and both are seen to be mutually consistent.Though estimates of λ u in the x direction exceed those in the y direction by about 25-30 %, we hesitate to interpret this as anisotropy due to their relatively large uncertainty.Figure 14a and b compare sample x-and y-directional variogram values, respectively, based on S 1 N , S 2 N and S 1 N and S 2 N jointly with variogram models calibrated against these values.
Whereas values of NLL corresponding to TPV models based on Gaussian and exponential modes are similar, those of KIC show a preference for exponential modes.Adopting the latter while considering S 1 N and S 2 N jointly yields λ u = 1.65 cm in the x direction and λ u = 1.31 cm in the y direction with an average of 1.48 cm.These correspond to ratios µ = λ u /L of upper cutoff to block size L equal to 0.055 in the x direction and 0.044 in the y direction with an average of 0.049.Corresponding estimates of the lower cutoff λ l are 6.8 × 10 −2 cm in the x direction and 1.2 × 10 −1 cm in the y direction with an average of 9.4 × 10 −2 cm.Adopting the suggestion of Di Federico and Neuman (1997) that µ = λ u /L = λ l / l m yields support (measurement) scales l m = 1.24 cm in the x direction and 2.69 cm in the y direction with an average of 1.88 cm.The latter is about 12 times the inner radius of the MSP.Albeit one should consider all the approximations involved in this estimate, we note that it is consistent with a definition of MPS support volume by Tartakovsky et al. (2000) as a region containing 90 % of total gas flow (see their Fig. 6).Estimates of µ and l m for all cases are listed in Table 2.
2. Identification of this nonlinear power-law scaling was greatly enhanced by a method of analysis that extend its range to virtually all lags (Guadagnini and Neuman, 2011) known as Extended Self-Similarity (ESS) and a generalized version thereof (G-ESS).
3. Our estimates of the Hurst scaling exponent were found to increase with support scale, implying a reduction in roughness (anti-persistence) of the log permeability field with measurement volume.
4. Tidwell and Wilson (1999) were able to characterize log permeabilities associated with all tip sizes by stationary variogram models.This, coupled with our findings that log permeability increments associated with the smallest tip size are approximately Gaussian and those associated with all tip sizes scale in the manner of multifractals at intermediate lags, are consistent with a view of the data as a sample from truncated fractional Brownian motion (tfBm).Since in theory the scaling exponents, ξ(q), of tfBm at intermediate lags vary linearly with q we conclude, in accord with Neuman (2010aNeuman ( ,b, 2011)), that nonlinear scaling in our case is not an indication of multifractality but an artifact of sampling from tfBm.

Our demonstration in Appendix
A that tfBm is consistent with ESS scaling according to Eq. ( 6) at all separation scales, and with power-law scaling according to Eq. ( 2) at intermediate scales, explains why and how ESS works for our data at all scales.The same explains how and why ESS worked for sub-Gaussian processes W G(s; λ l , λ u ) considered by Guadagnini and Neuman (2011).
6.The fact that our data are consistent with Eq. ( 6) but not with Eq. ( 2) at small and large lags constitutes yet another indication that, despite their nonlinear power-law scaling at intermediate lags, the data are inconsistent with multifractals or fractional Laplace motions, which theoretically scale in this manner at all lags.The same likely holds true for other Gaussian or heavy-tailed earth and environmental variables (such as those listed in our introduction) that scale according to Eq. ( 2) at intermediate lags and according to Eq. ( 3) over an extended range of lags, a possibility noted earlier by Guadagnini and Neuman (2011).
7. Since increments of Y associated with the smallest tip size have a near-Gaussian distribution, we were able to identify the functional form and estimate all parameters of the corresponding tfBm based on sample structure functions of first and second orders.Our estimate of lower cutoff is consistent with a theoretical support scale of the data.
1 A is a coefficient, H is a Hurst exponent (0 < H < 1), i = 1 for tfBm with modes (defined in Neuman, 2010a) having exponential autocorrelation functions and i = 2 for modes having Gaussian autocorrelation functions.Figure A1 compares TPVs based on Gaussian modes with A = 1, H = 0.3, λ l = 10 and four values of λ u ( 10 4 , 10 3 , 5 × 10 2 , 10 2 ) with a power variogram (PV) γ 2 (s) = A 2 s 2H where A 2 = A(π/4) H (1 − H )/2H .The slopes of the TPV and PV coincide in a midrange of lags (labeled Zone II) but not in the outlying ranges of small and large lags (labeled Zone I and III, respectively).This break in power-law scaling at small and large lags is due entirely to the presence of lower and upper cutoffs, respectively, being unrelated to noise or oversampling which play no role in Fig. A1 (Neuman, 2010a).It follows that estimating H as the slope of the variogram on log-log scale is valid at intermediate lags but not at small and large lags which would lead, respectively, to over-and underestimation of its value.Figure A2 complements this analysis by juxtaposing the TPVs associated with Gaussian modes in Fig. 14a with corresponding PVs.For a PV Eq. (A2) takes the form S q = | G(s; λ l , λ u )| q = (q − 1)!! 2 A i q s qH 2 π if q is odd 1 if q is even q = 1, 2, ... n (A8) rendering a log-log plot of S q versus s linear with constant slope qH.As in the case of q = 2, the slopes of corresponding truncated structure functions are similar in the midrange of lags but larger and smaller, respectively, at small and large lags.

Fig. 4 .
Fig.4.Sample structure functions of absolute increments of various orders q versus lag along x direction on face 1 and (a) r i = 0.15 cm, (b) r i = 0.31 cm, (c) r i = 0.63 cm, (d) r i = 1.27 cm.Dashed vertical lines delineate ranges of lags within which power-law scaling is noted.
Figure 5 Fig.5toward linear variation with q.In Appendix A we explain this behavior theoretically by demonstrating that tfBm scales according to Eq. (2) at intermediate lags and according to Eq. (3) at all lags.The fact that our data scale according to Eq. (2) at intermediate lags allows us to follow an approach patterned afterGuadagnini and Neuman (2011): adopt the value of ξ(1) from Fig.5as computed by the method of moments, fit straight lines by regression to log S q versus log Figure 7 Figure 8 Figure 9 Figure10

Fig. 11 .
Fig. 11.Sample structure functions of absolute increments of various orders q versus lag in z direction and (a) r i = 0.15 cm, (b) r i = 0.31 cm, (c) r i = 0.63 cm, (d) r i = 1.27 cm.

Fig. 14 .
Fig. 14.Variograms obtained from multiple faces data in (a) x (faces 1, 2, 4) and (b) y (faces 1, 3, 5) directions, on the basis of S 1N ( ) and S 2 N (×).Estimated variograms are also reported with continuous (Gaussian modes) and dashed (exponential modes) lines.Black lines correspond to estimated variograms obtained on the basis of S 1 N and S 2 N jointly.

Table 1 .
Calibration results, main statistics and Model quality criteria.The 95 % confidence intervals of the parameter estimates are reported in parenthesis.

Table 2 .
Multiple faces data.Estimates of µ = λ u /L and the associated support scale l m .