Multifractal detrended fluctuation analysis in examining scaling properties of the spatial patterns of soil water storage

Knowledge about the scaling properties of soil water storage is crucial in transferring locally measured fluctuations to larger scales and vice-versa. Studies based on remotely sensed data have shown that the variability in surface soil water has clear scaling properties (i.e., statistically self similar) over a wider range of spatial scales. However, the scaling property of soil water storage to a certain depth at a field scale is not well understood. The major challenges in scaling analysis for soil water are the presence of localized trends and nonstationarities in the spatial series. The objective of this study was to characterize scaling properties of soil water storage variability through multifractal detrended fluctuation analysis (MFDFA). A field experiment was conducted in a sub-humid climate at Alvena, Saskatchewan, Canada. A north-south transect of 624-m long was established on a rolling landscape. Soil water storage was monitored weekly between 2002 and 2005 at 104 locations along the transect. The spatial scaling property of the surface 0 to 40 cm depth was characterized using the MFDFA technique for six of the soil water content series (all gravimetrically determined) representing soil water storage after snowmelt, rainfall, and evapotranspiration. For the studied transect, scaling properties of soil water storage are different between drier periods and wet periods. It also appears that local controls such as site topography and texture (that dominantly control the pattern during wet states) results in multiscaling property. The nonlocal controls such as evapotranspiration results in the reduction of the degree of multiscaling and improvement in the simple scaling. Therefore, the scaling property of soil water storage is a function of both soil moisture status and the spatial extent considered.


Introduction
The spatial and temporal pattern of soil water storage is an important input variable in assessing land-atmosphere interactions, infiltration, recharge, and performance of engineered covers.Water storage is also a key input in monitoring the soil water balance and validation of several models (Rodriguez-Iturbe et al., 1999).The variability in soil water storage is shown to be strongly related to topographic, geologic, soil, and vegetation parameters (Braud et al., 1995;Moore et al., 1988).These physical factors and environmental processes (rainfall, evapotranspiration, runoff, and snow melt) do not operate independently, but as an ensemble of processes with a complex and nested effects.This, in turn, results in a pattern of soil water storage that varies as a function of spatial scale.Several studies have reported a scale dependent pattern and variability of soil water storage (Kachanoski and de Jong, 1988;Gomez-Plaza et al., 2000;Kim and Barros, 2002;Biswas and Si, 2011a).
In order to examine how information transfers from one scale (e.g.pedon scale) to another (e.g.satellite image scale); we need the scaling characteristics of soil properties.Scaling analyses such as fractal and multifractal analyses require the data series to be stationary.Hu et al. (1997Hu et al. ( , 1998) ) and Rodriguez-Iturbe et al. (1995) have characterized the apparent disorder in spatial organization of soil moisture and reported that the variance of soil moisture follows power law decay, typical of a scaling process.If such scaling laws are found to be appropriate in describing the fluctuation of soil moisture over wider range of scales, then the characterization of soil moisture could be based on the theory of statistical self-similarity that provide linkage over a range of spatial scales (Hu et al., 1998).However, these studies are based on spatial scales (i.e., resolutions) of 30 × 30 m 2 and only for the surface <5 cm soil depth.Such information can be utilized in specific land management practices and precision agriculture provided that the observed scaling law holds at a field scale and to a deeper soil depth.This will also open a window of opportunity during aggregation of point measurements and disaggregation of remotely sensed data (Lin, 2003).However, a priori assumption of similarity in scaling laws between the fluctuations in the remotely sensed large scale and small field scale soil water storage data could be erroneous owing to the increased importance of localized factors as the resolution increases and need separate treatment.
Spatial variability includes spatial trends (non-stationary) and fluctuations (stationary).Several studies have indicated that non stationarities introduce superficial scaling features (Bhattacharaya et al., 1983;Koscielny-Bunde et al., 2006).Such scaling features are not intrinsic to the variables under study because they are artifacts of erroneously treating nonstationary fields as stationary.Nor do they reflect the actual in situ fluctuation that is the result of numerous co-existing and interacting factors and processes, and hence, scale transformations based on such relationships could be erroneous.To remove the undue influence of the larger scale trends on the scaling properties at a scale of interest, there is a need to identify the intrinsic scaling property of the fluctuation in soil water storage at a field scale that results from interaction of all the underlying processes.Also there is a need to evaluate the presence as well as extent of spatial scales with a particular scaling property.Such analysis is useful in defining the spatial extent over which simple scaling up of point observations is possible as well as predicting soil water based on the observed scaling laws.
Soil water distribution within heterogeneous fields is often complex owing to the numerous physical factors and processes controlling its spatial and temporal variability.The large scale processes results in long-range correlations, with an autocorrelation (two-point correlation) function decaying slowly with increase in separation distance.If the two-point correlation function decays as a power-law, we have a scaling phenomenon.The power spectrum P (f ), which is the Fourier transformation of the correlation function, will also be a power-law function: P (f ) = f −c , where f is the frequency and c is the scaling exponent.However, the twopoint correlation function may not be the best way to characterize the long-range correlation, because it does not take into account of the structures in between two points, and it does not measure correlation between two units larger than a point.As the long-range correlation is very common in nature, engineering, and medicine, many methods have been developed to analyze the scaling property in the long-range correlation.Proven methods that remove nonstationarities from the data series include the Hurst rescaled-range analysis (Hurst, 1951), the wavelet transform modulus maxima (WTMM; Arneodo et al., 2002;Zeleke and Si, 2007), and the detrended fluctuation analysis (DFA) (Peng et al., 1992(Peng et al., , 1994)).The Hurst re-scaled range analysis is based only on the first moment of the series and hence does not provide the detailed characterization of the spatial statistics, nor does it remove nonstationarities (Liebovitch et al., 2002;Koscielny-Bunde et al., 2006).Similar to the wavelet transform (Mallat, 1999), the WTMM is easily applied to binary observations (i.e., n = 2 k , where n = number of observations and k = 1, 2, 3, 4...) for its computational simplicity and to minimize "edge effects" (Mallat and Hwang, 1992;Oświe ¸cimka et al., 2005).The DFA, on the other hand, is more computationally simple, straight forward and flexible technique and does not depend on any specific number of observations (e.g., binary observations) given it covers the scales of interest (Kantelhardt et al., 2002;Oświe ¸cimka et al., 2005;Koscielny-Bunde et al., 2006).Integration of the data series reduces outliers that usually exist in spatial data whereas the detrending procedure reduces the effect of nonstationarities.The DFA technique, originally developed for monofractal variables, has been extended to accommodate analysis of scaling heterogeneity through multifractal detrended fluctuation analysis (MFDFA; Kantelhardt et al., 2002).This development was achieved through modification of the traditional Hurst function, H into a generalized function, h(q) so that moments of different orders can be evaluated.A detailed comparison between WTMM and MFDFA can be found in Kantelhardt et al. (2002).
For improved modeling or prediction of soil water distribution we need to understand its response to these numerous factors and processes than just to a single extrinsic factor that introduces monotonous trend.The presence of such trend can obscure the scaling property of the series that reflects both local and global processes.It is also important to minimize the effect of measurement artifacts that introduce outliers or noise that could be easily propagated during upscaling of observations.The DFA is a noble tool that handles these problems so that final output reflects the effect of the suite of processes on soil water distribution.The DFA technique has been successfully applied for scaling studies of nonstationary series in climatology (Kurnaz, 2004;Kiraly and Janosi, 2005), river runoff (Koscielny- Bunde et al., 2006), financial series (Grau-Carles, 2001), and biomedicine (Peng et al., 1992).However, to the best of our knowledge, this noble technique has not been applied in characterization of scaling properties of soil water storage.Utilizing this noble technique, we will explore scaling properties of soil water storage under field conditions.The specific questions we will address: (i) Are there simple scaling (monofractal type) and multiscaling (multifractal type) characteristics in the spatial correlations of soil water storage?(ii) Is there any change in the scaling pattern with the change in climatic processes?

Scaling analysis based on the detrended fluctuation series
The DFA for a one dimensional data series can be described as follows.Let X j be a data series of length N with a compact support (j = 1, 2,•••, N).The integrated series or the profile at location i, Y (i), is determined by taking the sum of deviations from the mean value (Kantelhardt et al., 2002;Telesca et al., 2004) i.e., where X is the mean value of the series for N observations.The integrated series Y (i) is then divided into N s nonoverlapping segments of equal size s.Since the length N of the series may not be a multiple of s, an unequal and short part (< s) of the profile may left at the end.In order not to disregard this part of the series, the same procedure is repeated starting from the other end.Thus, 2N s segments are obtained altogether.Then the local trend for each of the 2N s segments is calculated by a least squares fit of a polynomial function.
The mean squared difference between the series [Y (i)] and the ordinate of the fitted polynomial [y v (i)] is calculated as Note that the indices i and v correspond to the original data points and the segments of size s, respectively.The F (s,v) is used to calculate the fluctuation functions as follows.The standard fluctuation function, F 2 (s) is calculated as a square root of the average of F (s,v) over all the segments of size s, i.e., The standard fluctuation function is based on the variance (second order moment) of the fluctuation of observed values relative to fitted polynomial trends.The fluctuation function can be extended to include higher order moments (say q values) to analyze the scaling property of different ranges of fluctuations, and also the detrending polynomial, y v can take any order n (linear, quadratic, cubic, etc.).The generalized fluctuation function, F q (s) is thus defined as where the variable q can take any real value except zero.In the case where q is zero, the fluctuation function cannot be determined directly from Eq. ( 4) because of the diverging exponent.Thus, F 0 (s) is approximated by taking the logarithmic average as, Repeating the above procedure for several length scales, a relationship can be developed between the fluctuation function and the segment length.Typically, F q (s) will increase with increase in scale s.If the series X i has a long range power law correlation, F q (s) increases with increase in s as a power law where h(q) is the generalized Hurst scaling function (Telesca et al., 2004).For q = 2 we have the standard DFA analysis.
In this case, the scaling exponent h(2) provides information about the average fluctuation of the series.The series can be categorized into one of the following three types depending on the h(2) value.These are: (i) 0 < h(2) < 0.5 for an anti-persistent type long range correlated process where large values (compared to the average) are more likely followed by small values and vice versa, (ii) h(2) = 0.5 for an entirely random uncorrelated distribution, and (iii) 0.5 < h(2) < 1 for a persistent and long range correlated process where large values are more likely to be followed by large values and vice versa.

Scaling heterogeneity (multifractality) of the detrended fluctuation series
The link between the generalized fluctuation function and the standard box counting formalism of multifractal analysis is also a straight forward one (Kantelhardt et al., 2002).For a normalized series X k , the mass distribution probability in the v-th segment of size s unit, P s (v), can be calculated as The mass scaling function, τ (q) is then defined via the partition function µ(q,s) as, The mass scaling function is related to the generalized Hurst scaling function, h(q) as (Kantelhardt et al., 2002) τ (q) = qh(q) − 1 .
A recent study by Yu et al. (2011)  therefore needs more in-depth and rigorous testing of the proposed relationship, where H is the non-conservation parameter in the universal multifractal formalism.Therefore, in this manuscript we have used the well established relationship between τ (q) and h(q) proposed by Kantelhardt et al. (2002).
The singularity spectrum of a multifractal measure, f (α), is related to τ (q) via a Legendre transform as (Feder, 1988) where, α is a singularity strength or Hölder exponent and obtained numerically as the first derivative of the τ (q) function with respect to q. Hence f (α)denotes the dimension of the subset of the series that is characterized by α herefore, f (α) can be used during the multifractal analysis for convenience and ease of interpretation once the series is normalized, integrated, and detrended using the appropriate detrending polynomial.
Schertzer and Lovejoy (1987) derived a universal multifractal (UM) model based in certain reasonable assumptions about the mechanism generating multifractals.The critical assumption was that the underlying generator is a random variable with an exponentiated extremal Lévy distribution (Zeleke and Si, 2006).The UM model can describe the τ (q) function under the assumption of conservation of mean value of the variable.
where, α (commonly known as Lévy index) indicates the degree of multifractality based on the deviation of τ (q) function from a monofractal type of scaling (UM model) (Seuront et al., 1999) and C1 indicates the co-dimension between observation space and fractal dimension.Nonstationarity is a common aspect of complex variability and is often be associated with different trends in the data series or patches with different local statistical properties (Kantelhardt et al., 2001).The DFA approach described above reduces the effect of such nonstationarities on scaling property of a variable.The reason for detrending analysis is to remove the undue influence of larger (than the scale of interest) scale on the statistics of soil water at the scale.At small scale, we remove trends of small scale; at large scale, we remove only large scale trend and small scale trend remains, because the small scale trends do not affect the scaling properties at a large scale.Therefore, we do not lose any critical information in soil water content.This is evident from the following three key aspects of the method.Firstly, the original series (i.e., as determined by all hydrological processes) is integrated into a continuous profile by summing up the deviations from the mean value (see Eqs. 1 and 2).In other words, the integrated profile of the series is the result of all the processes determining the spatial pattern of the variable.Secondly, the integrated profile was divided into segments and regression lines were fitted to each of these segments (window).When the fitted trends are removed, what remains is a fluctuation function, which is the difference between the integrated series and its best fit regression line at a given window size.This function contains all the local as well as global feature of the data series which is free from nonstationarities.Thirdly, the size of the segment (the window) was continuously varied later on; i.e., during the scaling analysis.At this step, it is important to note that the scaling property is determined by relating the fluctuation function to several window sizes (Eq.6).These points clearly show that the DFA procedure does not exclude any process that determine the hydrology of the system; rather transform resulting data series into a fluctuation function where the effect of nonstationarities is significantly reduced.Thus, in essence, the main advantage of the DFA method is that it allows detection of scaling property of a physical variable that is embedded in a noisy data or containing monotonous polynomial trends that can mask true fluctuations of the series.

Material and methods
The study site is located in a sub-humid climate at Alvena, Saskatchewan, Canada.The geographical location of Alvena is 49 • 44 N latitude and 107 • 35 W longitude.The site has rolling topography (locally referred to as a hummocky terrain) with a dominant soil type of an Aridic Ustoll (US Soil Taxonomy).The surface texture is loam to clay loam.Long term mean annual air temperature is 2.2 • C and precipitation is 350 mm.The potential evapotranspiration is 624 mm yr −1 , resulting often in a water deficit of more than 250 mm yr −1 .A north-south transect of 624 m length was established in 2001 and soil water storage have been monitored at 104 locations (at 6 m regular intervals) at several times (more than 30 series) using gravimetric and capacitance probes.For this particular study, six series of soil water storage data (between 2002 and 2005) of the surface 0 to 40 cm depth was selected.To reduce edge effect, only 93 locations were selected for the data analysis.Six water storage series were selected based on their representativeness of the dominant climatic processes of the area that determine soil water storage pattern (snow melt, rainfall, and evapotranspiration) and measurement methods used (only gravimetrically measured data was used for better accuracy).These were 14 August 2002, 14 September 2002, 23 May 2003, 23 October 2003, 1 August 2004, and 30 May 2005.This manuscript focuses on scaling properties for soil water storage patterns of the top 0 to 40 cm depth.Changing spatial scale, as used in here, refers to varying the number of discrete data points (of fixed spacing) used in deriving the mean value.This can also be referred to as "coarse graining" of discrete samples, and integrates the concept of both support and spacing described in Western and Blöschl (1999).For instance, taking the mean value of four consecutive discrete samples that are initially spaced 3 m apart is equivalent to increasing the scale from 3 m to 12 m.In this transform both the spacing and support are upscaled.In other words, we are increasing both the spacing and support (i.e., making a practical assumption that the mean value is representative of the length within the spacing).
Water content was measured following the standard gravimetric technique (Gardner, 1986).The core samples were collected from each measurement locations using a 5 cm (internal) diameter sampler that was vertically mounted on a truck with a hydraulic system.The undisturbed core samples were sliced into 10 cm increments and placed in a plastic bag for oven drying.The mass based water content is determined as mass difference between the wet and oven-dry samples divided by the sampling volume.The bulk density of the sample is determined as the oven dry mass of the soil sample divided by the sampling volume.The volumetric equivalent of the gravimetric water content of the samples was determined by multiplying the mass based water content with bulk density.Mean values of four core samples, taken at 10 cm vertical intervals, was used to obtain a data point representing the surface 0 to 40 cm depth.Soil water storage in the root zone is critical to plant productivity.The top 60 cm have majority of roots.However, soil water content at depth determined using gravimetric methods may not be accurate because of the potential compaction when a punch truck is used to extract soil cores in the clay soil.Our experience in this field suggested that 40 cm depth is free of observable compaction and therefore, we chose 40 cm.
Particle size distribution of the soil samples were determined from the first measurement sets using the hydrometer method (Gee and Bauder, 1986).Organic carbon content was determined using LECO-12 carbon determinator (LECO Corporation, St. Joseph, MI).Relative elevation along the transect was determined using a Laser Theodolite Total Station (Sokkisha Electronic Total Station, Set 5, Sokkisha, Tokyo, Japan).Both the DFA and MFDFA analysis were performed using programs written in Mathcad Professional (version 12, Mathsoft Inc., Cambridge, MA, 2002) and Statistical Analyses Software-SAS Version 8 (SAS Institute Inc., Cary, NC).

Water storage series and order of detrending polynomials
Figure 1 shows the spatial distribution of soil water storage at six occasions, two soil properties (clay and organic carbon) and one topographic variable (relative elevation).The monthly mean precipitation values for the years 2002 to 2004 are shown in Fig. 2. The total precipitation received in the calendar years of 2002 and 2004 was generally higher than the long term average precipitation of the area and these years are regarded as wet years.However, high soil water storage was observed only in May 2003 and 2005.The observed high soil water in the two series compared to others appears to be the result of relatively high snow fall in the months preceding the measurements and gradual melting of snow where downward infiltration significantly exceeds run off.Water storage during all the occasions were negatively correlated (r 2 = 0.18 to 0.40; significant at p = 0.01) to relative elevation (RE) as expected (Table 1).The relationship to clay content (CL) was significant (r 2 = 0.22, 0.15, and 0.14; significant at p = 0.01) only for measurements taken in August 2002, September 2003, and August 2004.The variance in May 2003 and 2005 were higher than the other four series and appear to be caused by non uniform snow redistribution and runoff related spatial variations.Although there were some significant (p = 0.01) positive relationships to CL, the distribution of water storage was dominantly controlled by RE.Consequently, there were systematic and localized trends in the spatial distribution of water storage at all scales that follows the landscape pattern.Prior to scaling property analysis, nonstationarities due to such localized trends have to be removed (or minimized) through transferring the series into its fluctuation function.To this end, identifying the correct order of the detrending polynomial is the initial step.
As shown in Fig. 1 both topography and clay have linear and higher order trends (trends that can be represented by linear lines and higher order polynomials).These trends are also reflected in soil water storage series.Therefore, the fluctuation functions after different detrending polynomials needs to be examined for remaining trends and noise in a data prior to scaling analyses.To this effect, different order of detrending were evaluated like linear (DFA1), quadratic (DFA2), cubic (DFA3), and fourth order (DFA4) polynomial.The coefficients of determination (for linear fits) of the double log plots of the fluctuation functions vs. segment sizes of successive detrending polynomial orders were compared.The comparison continued until the difference between the trends in the fluctuation functions of two successive orders is insignificant.The significance between different orders of polynomials was tested using "students t" test by comparing the means between each consecutive detrending functions.For our data, there were significant differences between the original data and the first order (linear) detrending function (based on Students t test result).There were no significant difference between first (linear) and second (quadratic) order of detrending, and hence the first order (DFA1) was selected as the best detrending polynomial function to obtain a stationary series for scaling analysis.In general, though the increasing order of the detrending polynomials enables us to remove considerable portion of existing nonstationarities, the significance between two orders of detrending help us to select the order.There is also a concern   a) 0.26 (a) 0.02 (a) Significant at p = 0.01, (b) Standard deviation about altering the intrinsic pattern of the series if the order of the detrending polynomial is too high (Kantelhardt et al., 2001).Bunde et al. (2002) reported that results are reliable only for certain orders, above which DFA yield the same type of behavior.Since the first order polynomial successfully removed majority of existing trends in our water storage data, it is reasonable to assume that the majority of the trends were of linear orders.Consequently, further scaling analysis was carried out using the first order detrended data series.

Evaluation of the DFA for scaling property
The standard fluctuation function (DFA2) was evaluated for power law relationships between fluctuations and scale.The fluctuations for all water storage series showed an almost exact power law increase with observation scales (Fig. 3).The coefficients of determination for a linear fit of the double-log plots of the series were between 0.99 and 1.00 (n = 21).Such power law relationships indicate the presence of scaling laws (Hu et al., 1997).

Multifractal analysis
The scaling analysis presented above is based only on the second order moment or the variance of the fluctuation function (i.e., q = 2).But in most physical and biological data the scaling property of low and high values (relative to the average) is often different.Such observations imply the need for multifractal analysis in which the scaling property is represented by an array of scaling exponents rather than by a single one.To this end, the scaling analysis has been extended by including higher and lower order moments (q values), i.e., in multifractal analysis (Eqs. 8,9,and 11).Mass exponents,τ (q) were derived from the fluctuation functions for q values between −20 and 20 and plotted against the q values (Fig. 4).A linear reference line (similar to monofractal type of scaling) (Fig. 4) was created following the UM model of Schertzer and Lovejoy (1987) to compare and characterize the observed scaling properties (Eq.12).A 1094.00 SS SSR = Sum or Squared difference of Residuals between the τ (q) of the data and the simulated monofractal type distribution using UM model, SS = a significant difference (p < 0.01) using Chi-square (χ 2 ) goodness of fit test between the τ (q) of the data and the simulated model data, df = degrees of freedom used for χ 2 statistics evaluation, and p = probability.
chi-square test for goodness of fit (at p = 0.01) indicate that all the soil water storage series are multifractal in nature as the mass exponent curve is quite different from the simulated monofractal type of scaling (Zeleke and Si, 2006).The sum of squared difference of residuals (SSR) of mass exponents and simulated monofractal type of scaling are summarized in Table 2.The SSR value for the soil water storage data series of 23 May 2003 and 30 May 2005 is way larger than the rest of the measurements indicating higher degree of multifractality compare to other measurements.The high precipitation during the year of 2002 and 2004 led to high soil water storage.The post snow melt period as controlled by the several local and non-local controlling factors (Grayson et al., 1997) affected the spatial distribution of soil water storage making it more heterogeneous in nature.With time, heavy demand of evapotranspiration by plant community reduces this heterogeneity and the degree of multifractality towards fall season.The slope of the τ (q) plots of water storage series were measured in two cases; a single (q = −20 to 20) and −20 to 20 at 1.0 increments).The solid line is a linear reference created following the UM model of Schertzer and Lovejoy (1987) passing through τ (q = 0).segmented (q = −20 to 0 and q = 0 to 20) and summarized in Table 3. Statistical significance of the difference between the variances under these two cases was evaluated using the F statistics (Press et al., 1992).The difference between the variances under these two cases (single and segmented) was significant (p = 0.01); implying that the τ (q) functions were significantly different from a linear function.A nonlinear τ (q) function means multiple scaling (Evertsz and Mandelbrot, 1992;Olsson and Niemczynowicz, 1996), which requires a hierarchy of scaling exponents (multiscaling) in order to accurately represent the scaling property.The degree of non-linearity of τ (q) function can give idea about the de-(q) In order to study the local scaling patterns, the multifractal spectrum [f (q) vs. α(q)] of six water storage series were calculated and presented in Fig. 5.The wider the spectrum (i.e., the higher the α max -α min value), the higher is the heterogeneity in local scaling indices of the variable and vice versa.The height of the spectrum, f (q), corresponds to the dimension (frequency distribution) of these scaling indices.Low f (q) values correspond to rare events (extreme values in the distribution), whereas the highest value of f (q) is the capacity dimension, which is obtained by assuming uniform distribution in all the segments.The spectra for the May 2002 have the widest range of α value (α max −α min = 1.05) indicating the most heterogeneous scaling indices or possibility of multiscaling.The spectra for the 30 May 2005 has also similar range of α value (α max −α min = 0.95) indicating multiscaling nature.The difference in microclimate, for example the difference in slope, concavities, soil texture, organic carbon content or the catchment area (Grayson et al., 1997) affected the distribution of water during snowmelt period resulted in the variability of soil water storage.The explanation of this scaling property requires numerous dimensions indicating multifractal nature of scaling.The high demand of evapotranspiration leading to a uniform drying process over time substantially reduces the variability of soil water storage pattern as indicated by the α max −α min value of 0.75, 0.65, and 0.55 for 14 August 2002, 1 August 2004 and 14 September 2002 respectively.The gradual decrease of α value over time indicates the reduction in the degree of multifractality.
The reduction in the variability or the increase in uniformity of soil water storage leads to the scaling property simple.The α max −α min value for 23 October 2003 is 0.70, which is a bit higher than the other series of similar time.This multiscaling nature of soil water storage might have existed from the higher precipitation during the year of 2002.
The scaling dimensions for 23 May 2003 and 30 May 2005 series vary from 0.65 to 1.70 and 0.60 to 1.55 respectively, which means that the representation of the scaling property of these variables requires numerous dimensions whose values are bound between 0.6 and 1.7; however that of other series requires 0.55 to 1.3, 0.75 to 1.4, 0.5 to 1.05, and 0.7 to 1.4 respectively for 14 August 2002, 1 August 2004, 14 September 2002, and 23 October 2003.The spectra of both May (2003 and2005) series has slight longer tail to the right of the maximum f (q) value, which is a characteristic of multifractal measure.Note that the right side of the spectrum corresponds to lower data values that are amplified by negative q values, and hence the right skewed feature is the result of more heterogeneity in the distribution of lower data values.
There are two sources of multifractality in time or spatial series as described in Kantelhardt et al. (2002).These are due to broad probability density distribution (long tailings) and differences in autocorrelation types.The multifractality observed in the water storage series appears to be the result of differences in the autocorrelation types for the small and large fluctuations.For the 23 May 2003 and 30 May 2005 series, the spatial variation in fractal dimensions is very high (Figs. 4 and 5) and, therefore, can be represented as multiple scaling pattern.The spatial variation in the fractal dimension gradually decreased over time.As discussed in the previous sections, this series is unique in that it is the result of a uniform drying process (evapotranspiration) and the variability (compared to the May series in the same year) was substantially reduced.Note that it is not possible to tell the differences between the May series with the other series based only on simple statistics such as mean and variance of the distribution.However, removal of nonstationarities and the subsequent scaling analysis showed the actual similarity and differences in terms of spatial scaling property.
The above results suggest that in any watersheds with significant topographic variations, the scaling property of soil water storage pattern may be different during dry and wet periods.The scaling property is determined by the number, type, and spatial extent of processes controlling soil water dynamics.During spring snow melt and after summer rainstorms, the spatial distribution of soil water is determined by several local and non local controls including slope, concavities, soil texture, organic carbon content, catchment area, and subsurface lateral flow (Grayson et al., 1997).Consequently, during these periods, the scaling property in soil water storage pattern becomes more heterogeneous (i.e., both as a function of scale and location within the landscape) resulting in a multifractal type distribution than during the drier periods.Based on a study using remotely sensed (large scale) soil water data in sub humid environment of Oklahoma, Kim and Barros (2002) also reported multifractality in soil water storage as a result of temporal evolution in wetting and drying regimes.The authors reported multifractal nature of soil water distribution at α-scale range (<10 km) as well as at β-scale range (>10 km), which exhibited multifractal to noise type scaling when the soil moisture levels are lower than field capacity (Kim and Barros, 2002).However, Mascaro et al. (2010) reported a multifractal scaling of soil water distribution at all domains in wet conditions using remotely sensed soil water measurement.The authors ascribed this variability to the signature of rainfall spatial variability (Mascaro et al., 2010).Generally the surface soil layer is exposed to various meteorological and environmental forcing such as rainfall, wind, solar radiation and become more dynamic than the deeper layers (Hu et al., 2010;Biswas and Si, 2011a, b).Moreover, the adequacy of plant roots in the surface also makes the surface soil layer dynamic.Therefore, the study of scaling properties in soil water content at the surface few centimeters (such as remote sensing measurement) is more complicated and is highly variable in nature (Biswas and Si, 2011a, b).On contrary, deep soil layers are less responsive to the changes in meteorological conditions (Hu et al., 2010), have less root activity (Cassel et al., 2000) and less disturbed soil structure (Guber et al., 2003;Pachepsky et al., 2005), which increase the buffering capacity of soil water changes in the deep layers and create an hydrological inertia (Martínez-Fernández and Ceballos, 2003) in soil water dynamics.Moreover, the rapid changes of soil water at the surface do not represent the actual changes in the vadosezone soil water storage.Therefore, it is difficult to differentiate actual wet and dry situations and the scaling property of soil water distribution at those situations.In this study we have considered soil water storage up to 40 cm, which is deep enough (field observation and experience) to exclude the highly dynamic nature and include majority of the root to understand the vadose-zone soil water dynamics.Generally plants take up more than 70 % of the water they need from the top 50 % of the root zone (Feddes, et al., 1978;Morris, 2006).Therefore, scaling properties of soil water storage up to 40 cm depth represented much more realistic situations.

Conclusions
We studied the scaling properties of the fluctuations in soil water storage in a sub humid climate of Saskatchewan using data series selected from a long term monitoring experiment.The selected series represent extreme soil water regimes (dry and wet) and also reflect the main hydrological processes in the region (snow melt, rainfall, and evapotranspiration).The data were analyzed using the multifractal detrended fluctuation analysis technique in order to characterize the intrinsic scaling property of soil water storage.The results showed a multiscaling property (multifractal type) over the entire scales for all soil water storage series.The degree of multifractality changes with the change in climatic processes.The highest scaling heterogeneity (multifractality) was observed for the series in May (i.e., after spring snowmelt or in wet period).This scaling heterogeneity gradually decreases over time showing a simpler scaling law towards the end of fall season (drier period).This multifractal scaling nature is mainly due to the heterogeneity in soil water storage pattern as affected by the micro climate during post snowmelt period.The high demand of evapotranspiration results in a uniform drying process which substantially reduces the soil water storage variability leading to a simpler scaling in nature.The implication is that the disaggregation of observations (e.g.remotely sensed large scale data to a field scale) for soil water storage based only on scaling laws could be erroneous during recharge periods, especially after spring snow melt.Therefore for adequate representation of the field scale variability, we need more sampling (monitoring locations) during wet periods than dry periods.

Fig. 4 .
Fig. 3 640 Fig. 3. Double -logarithmic plots of the standard fluctuation functions fitted to linear equation.In order to avoid overlapping of the plots (and hence difficulty in comparisons) a constant values have been added to the fluctuation functions of each series."Scale" refers to a period (cycles) in meters.

Table 1 .
Mean and standard deviation of soil water storage at the Alvena site during six observation dates and relationships to relative elevations (RE), organic carbon (OC), and clay content (CL).

Table 2 .
Sum or Squared difference of Residuals between the τ (q) of the data and the simulated monofractal type distribution using UM model for the soil water storage of six observation dates.

Table 3 .
Slope of the mass exponent functions τ (q) of the water storage series and their standard deviations.Significant difference between the variances during single line fit (reduced model) and the segmented fit (full model), df = degrees of freedom used for F -statistics evaluation, and p = probability.