Fractal behavior of soil water storage at multiple depths

Wenjun Ji, Mi Lin, Asim Biswas, Bing C. Si, Henry W. Chau, and Hamish P. Cresswell 2 1 Department of Natural Resource Sciences, McGill University, 21111 Lakeshore Road, Ste-Anne-de3 Bellevue, Quebec, Canada, H9X3V9 4 2 Department of Soil Science, University of Saskatchewan, Saskatchewan, Canada, S7N5A8 5 3 Department of Soil and Physical Sciences, Lincoln University, PO Box 85084, Lincoln, Christchurch, 6 New Zealand, 7647 7 4 CSIRO Land and Water, Canberra, ACT, Australia, 2601 8 * Correspondence to: A. Biswas (asim.biswas@mcgill.ca Phone: +1 514 398 7620; Fax: +1 514 398 9 7990) 10 11


Introduction
Knowledge on the spatial distribution of soil water over a range of spatial scales and time has important hydrologic applications including assessment of land-atmosphere interactions (Sivapalan, 1992), performance of various engineered covers, monitoring soil water balance, and validating various climatic and hydrological models (Rodriguez-Iturbe et al., 1995;Koster et al., 2004).However, high variability in soil is a major challenge in hydrology (Quinn, 2004) as the distribution of soil water in the landscape is controlled by various factors and processes operating at different intensities over a variety of extents (Entin et al., 2000).The individual and/or combined influence of these physical factors (e.g., topography, soil properties) and environmental processes (e.g., runoff, evapotranspiration, and snowmelt) gives rise to complex and nested effects, which in turn evolve a signature in the spatial organization (Western et al., 1999) or patterns in soil water as a function of spatial scale (Kachanoski and de Jong, 1988;Kim and Barros, 2002;Biswas and Si, 2011a).This complexity makes the management decision difficult at a scale other than that of measurement.Therefore, it is necessary to transfer variability information from one extent (e.g., pedon) to another (e.g., large catchment), which is called scaling.
The scaling of soil water is possible if the distribution of some statistical parameters (e.g., variance) remain similar at all studied scales.This feature, known as scale-invariance, means that the spatial feature in the distribution of soil water will not change if the length scales are multiplied by a common factor (Hu et al., 1997).Generally, the soil water will have a typical size or scale, a value around which individual measurements are centered.So the probability of measuring a particular value will vary inversely as a power of that value, which is known as the power-law decay, a typical principle of the scaling process.Now, as the spatial distribution of soil water follows the power-law decay (Hu et al., 1997;Kim and Barros, 2002;Mascaro et al., 2010), the spatial variability can be investigated and characterized quantitatively over a large range of measurement extents using the fractal theory (Mandelbrot, 1982).When the spatial distribution of soil water is the response of some linear processes, the scaling can be done using a single coefficient over multiple scales and the distribution shows monofractal behavior.However, the spatial distribution of soil water is the nonlinear response of multiple factors and processes acting over a variety of scales and therefore needs multiple scaling indices (multifractals) for quantifying spatial variability (Hu et al., 1997;Kim and Barros, 2002;Mascaro et al., 2010).
The multifractal behavior in the surface soil water as a result of temporal evolution of wetting and drying cycles has been reported from the sub-humid environment of Oklahoma by Kim and Barros (2002).Mascaro et al. (2010) reported the multifractal behavior of soil water, which was ascribed as a signature of the rainfall spatial variability.Though these measurements can provide a quick estimate of soil water over a large area, they are limited to very few centimeters of the soil profile.These studies reported the multifractal behavior of only the surface soil water indicating the superficial scaling properties.Surface soil layer is exposed to direct environmental forces and is the most dynamic in nature.The scaling properties of surface soil water can be used for land management practices, provided the observed scaling properties remain the same for the deep layers such as vadose zone or the whole soil profile.Understanding the overall hydrological dynamics in soil profiles requires information on the scaling properties and the nature of the spatial variability of soil water over a range of scales at deep layers as well (Biswas et al., 2012c).The information on the similarity in the nature of the spatial variability of soil water between the surface layer and deep layers may also help inferring about the soil profile hydrological dynamics.Therefore, the objectives of this study were to examine over time the scaling properties of sub-surface layers and their relationship with surface layers at different initial soil water conditions.We have examined the scaling properties of soil water storage at each layer and their trend with increasing depth from the surface (cumulative depth) over a 5-year period from a hummocky landscape from central Canada using the multifractal approach.The relationship between the scaling properties of the surface layer and the subsurface layers was also examined using the joint multifractal analysis.

Study site and data collection
A field experiment was carried out at St. Denis National Wildlife Area (52 • 12 N, 106 • 50 W; ∼ 549 m mean above sea level), which is located 40 km east of Saskatoon, Saskatchewan, Canada.The landscape of the study area is hummocky with a complex sequence of slopes (10 to 15 %) extending from differently sized rounded depressions to irregular complex knolls and knobs, a characteristic landscape of the North American prairie pothole region encompassing approximately 780 000 km 2 from north-central United States to south-central Canada (National Wetlands Working Group, 1997).Some of these potholes are seasonal in nature meaning to store water in the spring (wet period) and drying out during late summer and in fall season (dry period) (Fig. 1).Variable water distribution within the landscape and in different landform elements such as side slopes, knolls, and depressions support vegetation differently.For example, the large amount of stored water in depressions provide a luxurious supply of water to growing plants compared to knolls (Fig. 1).A transect of 128 points (576 m long) extending in the north-south direction covering multiple knoll-depression cycles was established in 2004 at the study site to examine the soil water variation at field scale.The sample points were selected at 4.5 m regular intervals along the transect to catch the systematic variability of soil water.Soil water measurements were carried out at every 20 cm depth down to 140 cm along the transect over the period of 2007 to 2011, among which the surface soil water (0 to 20 cm) was measured using vertically installed time domain reflectometry (TDR) probes and a metallic cable tester (model 1502B, Tektronix, Beaverton, OR), while deeper layers down to 140 cm were measured using a neutron probe (model CPN 501 DR Depthprobe, CPN International Inc., Martinez, CA) (Biswas et al., 2012a).Soil water content data were then multiplied by depth and added together to obtain the overall soil profile water storage so as to examine the fractal behavior of soil water storage (SWS) at different depths over time.A detailed description of the study site, development of the transect, measurement of soil water and the calibration of measurement instruments can be found in earlier publications from this project (e.g., Biswas et al., 2012a).

Data analysis
Various methods including geostatistics (Grego et al., 2006), spectral analysis (Kachanoski and de Jong, 1988), and wavelet analysis (Biswas and Si, 2011a, b) have been used to examine the scale-dependent spatial patterns of SWS.These methods generally deal with how the second moment of SWS changes with scales or frequencies.When the statistical distribution of SWS is normal, the second moment plus the average provide a complete description of the spatial series.However, for other distributions (e.g., left skewed distribution) higher-order moments are necessary for a complete description of the spatial series.For example, let us define the qth moment of a spatial series z as z q .In this situation, for a positive value of q, the qth moment magnifies the effect of larger numbers and diminish the effect of smaller numbers in z.While, on the other hand, for a negative value of q, the qth moment magnifies the effect of small numbers and diminish the effect of large numbers in the spatial series z.In this way, using variable moments, we can look at the effect of the magnitude of the data in a series and better characterize its spatial variability.

Statistical self-similarity or scale invariance
Soil water is highly variable in space and time.If the variability in the spatial/temporal distribution remains statistically similar at all studied scales, the SWS is assumed to be selfsimilar (Evertsz and Mandelbrot, 1992).Self-similarity, also called scale invariance, is closely associated with the transfer of information from one scale to another.We used the multifractal analysis to explore self-similarity or inherent differences in scaling properties of SWS in this study.

Multifractal analysis
On the spatial domain of the studied field, multifractal analysis was used to characterize the scaling property of SWS by statistically measuring the mass distribution (Zeleke and Si, 2004).The spatial domain or the data along the transect was successively divided into self-similar segments following the rule of the binomial multiplicative cascade (Evertsz and Mandelbrot, 1992).This method required that the two segments divided from a unit interval to be of equal length.With regards to a unit mass M (a normalized probability distribution of a variable or measured in a generalized case) relating to the unit interval, the weight was also partitioned into [h × M] and [(1 − h) × M], where h was a random variable (0 ≤ h ≤ 1) governed by a probability density function.Sequentially, the new subsets with their associated mass were equally divided into smaller parts.In this way, multifractal analysis was able to describe the scaling properties for the higher-order moments compared to semivariogram, which can only measure the scaling properties of the second moment.In a special case, if the scaling properties do not change with q, the spatial series can be identified as monofractal, when one scaling coefficient is enough to characterize scaling property of SWS.Generally, the multifractal analysis is good at measuring the highly fluctuated mass (box size) within a scale interval.This also provides physical insights at all scales regardless of any ad hoc parameterization or homogeneity assumptions in the analysis (Schertzer and Lovejoy, 1987).
For SWS spatial series, the scale-invariant mass exponent, was termed as τ (q) (Liu and Molz, 1997): where z was the SWS spatial series, x was the lag distance and the symbol ∝ indicated proportionality.The τ (q) is widely used in multifractal analysis.If the plot of τ (q) vs. q [or τ (q) curve] has a single slope (i.e., a linear line), then the series is a simple scaling (monofractal) type.If τ (q) curve is nonlinear and convex (facing downward), then the series is a multiscaling (multifractal) type.In this study, we used the universal multifractal (UM) model of Schertzer and Lovejoy (1987)  in mean value of SWS, this model simulated a cascade process with a scaling function in an empirical moment.It is thus used here to compare and characterize the observed scaling properties with a reference to the monofractal behavior.The goodness-of-fit between the τ (q) curves and the UM model was tested using the chi-square test.The sum of squared residuals (SSRs) between the τ (q) curve and the UM model was also calculated to test the deviation.The τ (q) curves over the range of q values (in this study −15 to 15 at 0.5 intervals) were fitted with a linear regression line (referred to as a single fit).The linear fitting of the τ (q) curves with q < 0 and q > 0 (referred to as segmented fit) was also completed.
The difference between the mean of slopes and segmented fits (for positive and negative q values) was checked using the Student's t test.
In a similar manner to Eq. ( 1), the qth-order normalized probability measure of SWS, µ(q, ε) (also known as the partition function), is proven to vary with the scale size: where ε is scale size in the ith segment and p i (ε) is the probability of a measure p i (ε) and measures the concentration of a variable of interest (e.g., SWS) by dividing the value of the variable in the segment to the whole support length (e.g., to the whole transect of length L units) (Meneveau et al., 1990;Evertsz and Mandelbrot, 1992).The mass exponent τ (q) was related to the probability of mass distribution of SWS.Moreover, the fractal dimension of the subsets of segments in scale size ε was measured by the multifractal spectrum f (q).When a coarse Hölder exponent (local scaling indices) of α was at the limit as ε → 0, f (q) was calculated as follows (Evertsz and Mandelbrot, 1992): and the local scaling indices, α, were given by Noting that f (α) was determined through the Legendre transform of the τ (q) curve: and Jensen, 1989).
The multifractal spectrum is a powerful tool in portraying the similarity and/or differences between the scaling properties of the measures (e.g., SWS).The width of the spectrum (α max -α min ) was used to examine the heterogeneity in the local scaling indices.The wider the spectrum, the higher the heterogeneity in the distribution of SWS and vice versa.Similarly, the height of the spectrum corresponded to the dimension of the scaling indices.The small f (q) values indicated rare events (extreme values in the distribution), whereas the largest value was the capacity dimension (D 0 ) obtained at q = 0.
In addition to the multifractal spectrum, [f (q) vs. α(q)], for many practical applications, we used models to incorporate a few selected indicators to describe the scaling property and variability of a process.One of the widely used models for multifractal measure was the generalized dimension.The generalized dimension was calculated as follows: when q = 1, D 1 was referred to as the information dimension (also known as entropy dimension), which provided information about the degree of heterogeneity in the measure distribution in analogy to the entropy of an open system in thermodynamics (Voss, 1988).If the value of D 1 is close to unity, it indicated the evenness of measures over the sets of cell size, whereas the value approaching 0 indicated a subset of scale in which the irregularities were concentrated.The D 2 , known as the correlation dimension, was associated with the correlation function and measured the average distribution density of the SWS (Grassberger and Procaccia, 1983).For a monofractal distribution, D 1 and D 2 tend to be equal to D 0 .A same value of D 0 , D 1 , and D 2 indicates that the distribution exhibits perfect self-similarity and is homogeneous in nature.Contrarily, in multifractal type scaling, the D 1 and D 2 tend to be smaller than D 0 , showing D 0 > D 1 > D 2 .Accordingly, the D 1 /D 0 value can be used to describe the heterogeneity in the distribution (Montero, 2005).When this value equals to 1, it indicated exact monoscaling of the distribution.

Joint multifractal analysis
While the multifractal analysis characterized the distribution of a SWS spatial series along its geometric support, the joint multifractal analysis was used to characterize the joint distribution of two SWS spatial series along a common geometric support.As an extension of the multifractal analysis, the length of the data sets was also divided into several segments of size ε.Two variables (P i (ε) and R i (ε) representing two spatial series of SWS) were used here to measure the probability of the measure in the ith segment, when Among them, α and β were the local singularity strength, which respectively represented the mean local exponents of P i (ε) and R i (ε) in the corresponding expressions above.The partition function for the joint distribution of P i (ε) and R i (ε) was calculated as follows (Chhabra and Jensen, 1989;Meneveau et al., 1990;Zeleke and Si, 2004): where the normalized µ was the partition function, q and t were the real numbers for weighting, and the aforementioned local singularity strength (coarse Hölder exponents) α and β were the function to q and t as well: To indicate the dimension of the joint distribution, the multifractal spectra f (α, β) were given by In fact, the joint partition function in Eq. ( 6) can be simplified to Eq. ( 2) when q or t is equal to 0. In this case, the joint multifractal spectrum was transformed to the multifractal spectrum with a single measure.When both q and t were 0, f (α, β) reached maximum and indicated box dimension of the geometric support of the measures.Pair value of α and β fluctuates with the change of variable q and t.Therefore, it is possible to examine the distribution of high or low values (different intensity levels) of one variable with respect to another by varying the values of q or t.As the joint multifractal spectra f (α, β) represent the frequency of the occurrence of certain values of α and β, high values of f (α, β) represent strong association between the values of α and β.The Pearson correlation coefficient was used to quantitatively describe their relations across similar moment orders.In addition, correlation coefficients between the surface layer and subsurface layers were used as well to examine the similarity in the scaling properties.Additionally, a contour plot was used to represent the joint distribution of a pair of variables by permuting similar values (highs vs. highs or lows vs. lows) of q and t.The bottom left part of the contour graph presents the joint distribution of high data values of both variables while the top right part represents the low data values of both variables.Therefore, a diagonal contour with low stretch indicates a strong association between the variables in consideration (Biswas et al., 2012b).

Spatial pattern of soil water storage at different depths
Average SWS for the surface 0-20 cm layer over the 5year period was 5.51 cm.A slight decrease in SWS was observed at the immediate deep layer (20-40 cm) and a gradual increase thereafter.The 5-year average SWS was 5. 45, 5.48, 5.56, 5.61, 5.69, and 5.77 1).These top and bottom boundaries formed a wider range (3.76 cm) of the average SWS at the surface layer compared to that at the deepest layer (1.27 cm).A big range (2.00 cm) in the standard deviation (maximum = 2.43 cm and minimum = 0.43 cm) of the measurement at the surface layer (0-20 cm) was also observed compared to that at the deepest layer (120-140 cm; maximum = 1.28 and minimum = 0.76).This indicated large variations in SWS at the surface layer that gradually decreased at deeper layers.The coefficients of variation (CVs) at the surface layer (0-20 cm) varied from 10 to 43 % and at the deepest layer (120-140 cm) varied from 13 to 23 % (Table S1 in the Supplement).
The maximum SWS at the surface layer also varied widely (maximum = 13.96cm and minimum = 4.64 cm) compared to the deepest layer (maximum = 9.81 cm and minimum = 6.71 cm) (Table 1).There was a gradual decrease in the maximum value and increase in the minimum value from the surface to the deepest layer.The maximum SWS at different layers was more localized.For example, there was high SWS at different layers at the locations of 100 to 140 m and 225 to 250 m from the origin of the transect.These locations had very high SWS compared to the field-average because they were situated in the depressions while low SWS was observed on the knolls.
The variations in SWS with time were evaluated within a year.There was little change in the average SWS over measurements within the years from 2007 to 2011 except for 2008 (Table 1).For example, average SWS was 6.47, 6.03, 6.54, and 6.33 cm on 6 April 2010, 19 May 2010, 14 June 2010, and 28 September 2010, respectively.However, the average SWS in 2008 drops from 6.28 cm on 2 May 2008 to 3.51 cm on 17 September 2008 in the surface 0-20 cm layer.This falling trend was observed at all soil layers.When compared between years, the trend over time and with depth was very similar in 2007 and 2009 while slightly different between 2010 and 2011 (Table 1).A decreasing trend of the variability was also observed with time.For example, the CV of the surface layer was around 28 % on 2 May 2008, which gradually decreased to around 13 % on 17 September 2008 (Table S1).The average water storage for soil layers with increasing depth was also calculated by adding the individual layers together.The time-averaged values of SWS were 10. 96, 16.44, 22.00, 27.61, 33.30, and 39.07 cm for the 0-40, 0-60, 0-80, 0-100, 0-120, and 0-140 cm, respectively (Table S2).The CV of the 0-20 cm layer was the highest during the wet period and gradually declined to the smallest during the dry period (Table S3).The variability also gradually decreased with depth.

Statistical scale invariance
The power-law relationships and the statistical scale invariance were evaluated using a log-log plot of the aggregated variance of SWS spatial series at different depths of soil layers and the level of disaggregation (or scales) at different q values or statistical moments.The linear relationship of the logarithm of the variance with scale indicated the presence of statistical scale invariance (Fig. 2).The scale invariance was observed for all measurements and at all depths, although only all depths of three selected dates were presented as an example.The coefficient of determination (r 2 ) for a Figure 3. Mass exponents for soil water storage spatial series measured at selected 20 cm soil layer down to 140 cm in 2008 for a range of q (−15 to 15 at 0.5 increments).The solid line is a linear reference created following the UM model of Schertzer and Lovejoy (1987) passing through (q = 0).linear fit (n = 7) was between 0.99 and 1.00 (significant at P = 0.001) for any measurement days and depths.A similar trend in scale invariance was also observed for SWS with increasing depth.

Multifractal analysis
The τ (q) curves for the surface layer displayed deviation from the UM model during the wet period (Fig. 3).A high SSR value was observed between the τ (q) curves and the UM model.Nonlinearity in the τ (q) curve was observed and the slopes of the segmented fit of the τ (q) curves were significantly different from each other.For example, the SSR values between the τ (q) curve and the UM model were 27.74 and 50.49for the surface layer (0-20 cm) on 2 and 31 May 2008, respectively.The slopes of the τ (q) curve for single fit were 0.97 and 0.96, respectively, for the surface layer of 2 and 31 May 2008 (Fig. 3).The slopes of the segmented fit for these measurements were 1.04 (q < 0) and 0.87 (q > 0) and, 1.06 (q < 0) and 0.82 (q > 0), respectively (Fig. 3; Table S4).
With the maximum deviation at the surface layer, the τ (q) curves gradually became very similar to the UM model with depth.The SSR value decreased considerably in deep layers.The slopes of the τ (q) curve (single fit) became almost at unity with no significant difference with the UM model.There was no significant difference between the slopes of the segmented fit.For example, the SSR value was 6.17, 4. 98, 8.80, 8.50, 8.86, and 6.16 respectively for the 20-40, 40-60, 60-80, 80-100, 100-120, and 120-140 cm layer of 2 May 2008.The slopes (single fit) for these layers were 0.99, 1.00, 1.01, 1.01, 1.00, and 0.99, respectively (Fig. 3).The slopes of the segmented fit were also very close to unity with no significant difference between them.
The SSR values gradually decreased and the slopes became almost at unity with increasing depth (Fig. 4).For example, the SSR values were 14.11, 9.31, 7.71, 6.86, 6.71 and 6.30 and the slopes (single fit) were 0.98, 0.99, 0.99, 1.00, 1.00, and 1.00, respectively, for 0-40, 0-60, 0-80, 0-100, 0-120, and 0-140 cm layer (Table S5).The slopes of the segmented fit for the τ (q) curve became almost the same as soil layers went deeper (Fig. 4).The linearity of the τ (q) curves was gradually strengthened and the SSR value gradually fell with the depth increase of soil layers at any time.A significant difference was observed between the slopes of the τ (q) curves in segmented fitting at the surface layer of the first three measurements in 2007 (Fig. S1 in the Supplement), two measurements in 2008 (Fig. 4), three measurements in 2009, and all measurements in 2010 and 2011 (Fig. S2).
A decreasing trend in the SSR value was also observed over time within a year.During the dry period, the slopes (single fit and segmented fit) became almost at unity with no significant difference (Table S6).For example, the SSR value was 14.12, 8.25, 1.30, 1.46, and 0.52 and the slope was 0.99, 0.99, 1.00, 1.00, and 1.00, respectively, for the surface layer (0-20 cm) of 21 June 2008, 16 July 2008, 23 August 2008, 17 September 2008and 22 October 2008 (Fig. 3).Similarly, a small SSR value and consistent slope were also observed at the deepest layer (120-140 cm).The SSR values of the 120-140 cm were 2. 47, 2.47, 3.31, 3.44, and 4.57, respectively, for the measurements on 21 June 2008, 16 July 2008, 23 August 2008, 17 September 2008, and 22 October 2008  (Table S6).The slope (single fit) for all these measurements was equal to 1.01 (Fig. 3).There was very little difference in the slopes of the segmented fits.
A significant difference in the slopes of the segmented fit was observed for the surface layer (0-20 cm) of three measurements in 2007 (17 July, 7 August, and 1 September; Fig. S1), and three measurements in 2009 (21 April, 7 May, and 27 May) (Table S4; Fig. S2).The difference became nonsignificant with depth and during other measurement times.The trend in deep layers over time was very similar to that of 2008.However, the trend in the SSR values and the slopes with time was different in 2010 and 2011 (Table S6).There was very little difference in the SSR values at different times of the year.For example, the SSR value for the surface layer (0-20 cm) was 20. 79, 27.18, 24.63, and 26.66 and the slope (single fit) was 0.97, 0.97, 0.97, and 0.97, respectively, for the measurements on 6 April 2010, 19 May 2010, 14 June 2010, and 28 September 2010 (Fig. 3).The slope of the segmented fit of the surface layer (0-20 cm) was significant for all measurements in 2010 and 2011.However, the trend with depth was similar to other years.
The height of the multifractal spectrum at different depths of measurement was very similar over time.The width of the spectrum (α max -α min ) varied with depth and time (Fig. 5).Generally, a comparatively large value of α max -α min was observed at the surface layer during the wet period and the value gradually became smaller with depth.For example, the value of α max -α min for the surface soil layer (0-20 cm) was 0.23 and 0.31, respectively, for the measurements of 2 and 31 May 2008 (Fig. 5).Meanwhile, the value of α maxα min for the soil layers of 20-140 cm with 20 cm increment was 0.15, 0.14, 0.19, 0.20, 0.20, and 0.18 for 2 May 2008 and 0.25, 0.19, 0.11, 0.14, 0.12, and 0.11 for 31 May 2008, respectively (Fig. 6).In the later part of the year, the width of the spectrum gradually decreased (Table S8).For example, the α max -α min values were 0.19, 0.16, 0.07, 0.08, and 0.05, respectively for the surface layer on 21 June 2008, 16 July 2008, 23 August 2008, 17 September 2008and 22 October 2008.Similar trend in values of α max -α min was also observed at deep layers (Fig. 6).
The trend of the α max -α min values in 2007 and 2009 was very similar to that of 2008 (Table S8).A higher value of α max -α min was observed in the first three measurements of 2007 (Fig. S5) and three measurements of 2009 (Fig. S6).However, the values in the surface layer (0-20 cm) in 2010 and 2011 were always higher compared to the deep layers (Fig. 6).There was no decreasing trend in values for the surface layer over time.For example, the α max -α min value was 0.21, 0.24, 0.21, and 0.22, respectively, for the measurements on 6 April 2010, 19 May 2010, 14 June 2010, and 28 September 2010 (Fig. 6).However, the trend in the α maxα min value of deep layers was similar to that of other years.A similar trend was observed for cumulative SWS with increasing depth over the years (Fig. 7).Generally, the value of α max -α min was also small with the highest in the 0-20 soil layers and gradually decreased with depth (Fig. 7; Table S9).
A very similar height of the f (q) curve for all depths and all periods indicated a consistent frequency distribution of the scaling indices (Figs. 6 and 7).Additionally, the position and the symmetry of the curve revealed the distribution of scaling exponents.A symmetric f (q) curve indicated uniform distribution of the scaling exponents.The left side of the spectrum corresponded to the large SWS that were amplified by the positive values of q, whereas the right side indicated smaller SWS that were amplified by negative q values.Symmetry leaning towards the left side during the early spring and in the surface layers in 2008 clearly showed the wider distribution of scaling indices and multifractal nature of the SWS (Fig. 6).While the shifting of the symmetry towards the right side clearly indicated less variable scaling indices and thus reduction of multifractal behavior.During the wet years of 2010 and 2011, the symmetry towards the left side indicated the variability in the scaling indices.This also persisted with depth.A similar trend was observed for different years at all depth layers (Fig. 7).
Generally, the D 1 and D 2 values for different depths of different measurements were very close to 1 (Fig. 8 and Table S10).In general, the D 1 value of the surface layers gradually increased with depth.Similarly, at any depth, the D 1 values gradually increased from the spring to the fall season through summer (Fig. 8).The highest variation in D values with q was observed in the surface layer and in the spring season and gradually decreased with depth and later part of the growing season.For example, the first three measurements in 2007 and 2009 presented high D values at high q values (Figs.S9 and S10).This high D value gradually decreased in the dry period of the year.For example, the D value with positive q was high in the surface layer of 2 and 31 May 2008 (Fig. 9), whereas it gradually decreased at the later part of the year (e.g., 17 September 2008).The trend with time and depth in 2007 and 2009 was very similar to that of 2008 (Tables S10 and S11).A consistently high D value was observed in the surface layer for all 2010 and 2011 measurements (Fig. 9).The trend in D values with depth in 2010 and 2011 was also similar to other years.A high value of D 1 and D 2 were also observed at all depth layers for all measurements (Fig. 10; Table S11).

Joint multifractal analysis
There were strong correlations between the scaling property of the joint distribution of the surface soil layer and the deep soil layers.The narrow width and the diagonally oriented contours between SWS measured on 22 October 2008 at 0-20 and 20-40 cm layers clearly demonstrate strong association between those two layers (Fig. 11).The correlation between the surface 0-20 cm and the deep layers on 2 May 2008 (wet period) was larger than 0.9 (significant at P = 0.001; Table 2).The highest correlation was observed  between those layers closest to each other.The correlations gradually increased over time and showed high consistency between different layers on 17 September 2008 (Table 2).A very similar trend was observed in other years.

Discussion
The amount of water stored in the soil is the result of the dominant underlying hydrological processes.Located in semi-arid climate, the study area receives about 30 % of the long-term annual average precipitation as snowfall during   9. Generalized dimension spectra of soil water storage spatial series measured at each 20 cm soil layer down to 140 cm in 2008 for a range of q (−15 to 15 at 0.5 increments).
winter months (Pomeroy et al., 2007).Generally, the depressions receive snow from surrounding uplands or knolls as redistributed by strong prairie wind (Pomeroy and Gray, 1995;Fang and Pomeroy, 2009).The snow melts within a short period of time during the early spring and contributes a large amount of water.The frozen ground restricts infiltration and redistributes excess water within the landscape with greater accumulation in depressions (Fig. 1) (Gray et al., 1985).Apart from the snowmelt, the spring rainfall also contributes to the water inflow in the landscape (Fig. 1).This created a spatial pattern of SWS that was almost a mirror image of the spatial distribution of relative elevation (Biswas and Si, 2011a, c;Biswas et al., 2012a).
In the spring, the sources of water loss were the deep drainage and the evaporation.As the loss of water through deep drainage in the study area was as low as 2 to 40 mm per year, occurring mainly through the fractures and preferential flow paths (Hayashi et al., 1998;van der Kamp et al., 2003), the major loss occurred mainly through evaporation from the surface of the bare ground and standing water in depressions.These processes lose a very small amount of water compared to the input of water in spring and early summer leaving the soil wet.Moreover, the surface soil with high organic matter content and low bulk density stored a larger amount of water than the deep layers where the organic matter gradually decreased and the bulk density increased.Reflecting the long-term history of vegetation growth in the landscape, the variability of organic matter content (CV = 41 %) may be one of the main factors of the high variability in surface layer SWS (Biswas and Si, 2011b).
As the vegetation developed in summer, strong evapotranspiration resulted in the lowest average SWS.High amounts of water in the depressions allowed grasses to grow faster and transpire more water compared to the knolls (Fig. 1).
Figure 10.Generalized dimension spectra of soil water storage spatial series from surface to different soil layers (cumulative storage) at 20 cm increment down to 140 cm in 2008 for a range of q (−15 to 15 at 0.5 increments).
For example, the aquatic vegetation growth within the depressions was as high as 2 m, while the grasses on the knolls grew to a maximum of 1 m tall.The uneven growth of vegetation and the high evapotranspirative demand in summer narrowed the range of SWS.In the soil where water is more available, evapotranspiration will be stronger while the less evapotranspirative demand will be shown in the relatively dry soil.As a result, the excessive water in the relatively wet soil will be offset by evapotranspiration, reducing the disparities between maximum and minimum values.This variable water uptake was visible in the growth of vegetation in the later part of the growing season as well (Fig. 1).The reduction in the range of SWS was the largest in the surface layer and gradually decreased at deeper layers.This is because the surface layer was exposed to various environmental forces.For example, plants can take up more than 70 % of the water they need from the top 50 % of the root zone (Feddes et al., 1978).This dynamic behavior of the surface layer exhausted readily available water and finally reduced the range in water stor- age.This decrease in range also happened in the later part of the growing season.
The multifractal and joint multifractal analyses explained the scaling behavior of SWS at different depths over time.The linearity in the log-log plot between the aggregated variance in SWS and the scale at all soil layers over time indicated that SWS behaved under scaling laws (Fig. 2).The near unity slope of the τ (q) curves and the insignificant difference from the UM model indicated a monofractal type scaling at all layers except the surface layer during the wet period (until mid-to late June) where a multifractal behavior led to a slight convex downward curve (Fig. 3).This was also supported by a significant difference between the slope of single and segmented fit in the surface layer during the wet period.
Generally during the wet period, excess water fills and drains macropores quickly and creates variations in SWS.Variations in the evaporation due to uneven solar incidence over micro-topography also triggered SWS variability in the surface layer.Additionally, the snowmelt and the release of water controlled by local (e.g., soil texture) and non-local (e.g., topography) factors also affected the spatial distribution of SWS, making it more heterogeneous in the wet period (Grayson et al., 1997;Biswas and Si, 2012).To the contrary, as depth increased, less impact of environmental factors tended to create less variability in SWS and exhibited a monofractal behavior, which was consistent with the uniform slope shown in Fig. 3.During the dry period or later part of the growing season, the SWS storage variability at all depths was small and exhibited monofractal behavior (Fig. 3).Accordingly, the deeper layers in the wet period and all layers in the dry period can be accurately represented by only one scaling exponent while the surface layer in the wet period may require a hierarchy of exponents.A similar trend was observed in SWS of cumulative depth layers (Fig. 4).Resulting from increasing buffering capacity of the deeper soil layers, the variability of cumulative SWS overlaid the multifractal W. Ji et al.: Fractal behavior of soil water storage at multiple depths nature of the surface layer, and finally exhibited monofractal behavior in general.
The scaling patterns of SWS at different depths and periods were further examined using multifractal spectrum [f (q) vs. α(q)] (Figs. 6 and 7).The degree of convexity was used to characterize the heterogeneity of scaling exponents or the degree of multifractality.Large values of α max -α min indicated stronger heterogeneity in the local scaling indices of SWS or cumulative SWS and vice versa.The largest value for the surface layer(s) in the wet period indicated the most multifractal behavior of SWS.However, the value decreased with depth and gradually converged in deep layers (Fig. 6).This decline manifested a conformity in the scaling behavior of SWS at deeper layers.Over time, the α max -α min value of the surface soil layer decreased and became very similar to that of deep layers.This indicated a reduction in the degree of multifractality for surface soil layers from the wet period to the dry period.A consistent α max -α min value for all depths during the dry period suggested the homogeneity and least multifractal nature of SWS.A similar behavior was observed in the cumulative SWS (Fig. 7).
To sum up, both the unity slope of the τ (q) curves (Figs. 3  and 4) and the degree of convexity of the f (q) spectrum (Figs. 6 and 7) jointly demonstrated that dynamic behavior of surface soil layers in the wet period made SWS highly variable and exhibited a multifractal nature, while less environmental forces and increased buffering capacity of deep layers led to a monofractal nature.As a result, multiple scaling exponents were required to characterize the variability of SWS in the surface layer during the wet period, while a fewer number of exponents were necessary for deeper layers during wet period or all layers during dry period.
The height of the spectrum, f (q) revealed the dimension or frequency distribution of the scaling indices (Caniego et al., 2003).A low height of f (q) curve indicated rare events or extreme values in the distribution, while a high value represented uniform distribution in all segments.A very similar height of the f (q) curve for all depths and all periods indicated a consistent frequency distribution of the scaling indices.
The two upper soil layers during the wet period tended to exhibit a longer tail of the curve on the left, showing more heterogeneity in the distribution of large values.However, when stepping into the dry period, the spectrum tended to display a longer tail on the right compared to the left side, suggesting more heterogeneity in the distribution of smaller values.A few locations with standing water lead to the spatial differences during the wet period, whereas a few points with very small SWS due to high evapotranspiration by growing vegetation during the dry period results in the heterogenic distribution in smaller values.
The generalized dimension, D q was subsequently used to characterize the scaling property and variability in SWS (Figs. 9 and 10).The largest value of f (q), referred to as the capacity dimension (D 0 ) obtained at q = 0, was close to unity for all layers at different times (Fig. 9).The information dimension (D 1 ) obtained at q = 1 was different from the correlation dimension (D 2 ), which is denoted as the average distribution density of the measurement for the surface layers in the wet period (Grassberger and Procaccia, 1983).In this case, the different values of D 0 , D 1 , and D 2 indicated the multifractal nature of the distribution of SWS.Similarly, a non-unity value of D 1 /D 0 (Montero, 2005) also indicated the multifractal nature of SWS at the surface layer(s) during the wet period.However, over the growing season, the D 1 and D 2 value approached to D 0 and indicated a monofractal type behavior.Similar values of D 0 , D 1 , and D 2 during the dry period also indicated homogeneous distributions.
Joint multifractal distribution between the surface to various subsurface layers indicated the similarity in the scaling patterns (Table 2).Basically, the hydrological processes of shallower layers were similar to those of the top layer, while deeper layers showed more disparities from the surface.The nearest subsurface (20-40 cm) layer showed generally the highest similarity with the surface (0-20 cm) layer.However, in the wet period, the subsurface layers displayed the smallest similarity to the surface layer, suggesting a higher dynamic nature of hydrological processes.In the dry period, a stronger effect of vegetation overwhelmed the effect of small variations of water distribution, thus creating a more uniform distribution of SWS at all soil layers (Table 2).
Overall, our result revealed a multifractal behavior of surface soil layers during the wet period due to the dynamic nature of hydrological processes.This behavior gradually changed with depth and time (Fig. 12).In the deeper layers during the wet period, the behavior became less multifractal or nearly monofractal.Similarly, in the dry period, the vegetation development and its high evapotranspirative demand in the semi-arid climate of the study area increasingly buffered the variation of SWS, as a result, all the soil layers showed uniform distribution or monofractal behavior (Fig. 12).

Summary and conclusions
The transformation of information on soil water variability from one scale to another requires knowledge on the scaling behavior and the quantification of scaling indices.Surface soil water can be easily measured (e.g., remote sensing) and presents multi-scaling behavior (requiring multiple scaling indices).However, land-management practices require the understanding of the hydrological dynamics in the root zone and/or the whole soil profile.
In this manuscript, the scaling properties of soil water storage at different soil layers measured over a 5-year period were examined using multifractal and joint multifractal analysis.The scaling properties of soil water storage mainly suggested a monofractal scaling behavior.However, the surface layer in the wet period or with high soil water storage tended to be multifractal, which gradually became monofractal with depth.With the decrease in soil water storage, the scaling behavior became monofractal during the growing season.In a year with high annual precipitation, the soil stored more water in the surface layer throughout the growing period and displayed nearly multifractal scaling behavior.This multifractal nature indicated that the transformation of information from one scale to another at the surface layer during the wet period requires multiple scaling indices.On the contrary, the transformation requires a single scaling index during the dry period for the whole soil profile.The scaling properties of the surface layer were highly correlated with those of the deep layers, which indicated a highly similar scaling behavior in the soil profile.The study was conducted in an undulating landscape from a semi-arid climate and the results were very consistent over the years.Therefore, the observation completed at the field scale in this type of landscape and climate may be generalized in similar landscapes and climatic situations, otherwise may need to be examined thoroughly.The method used here can be transferred to examine the scaling properties in other experimental situations.
The Supplement related to this article is available online at doi:10.5194/npg-23-269-2016-supplement.

Figure 1 .
Figure 1.Conceptual schematics showing the vegetation growth patterns over the landscape at different times of the year.The figure is developed based on field observations and the scale is arbitrary.

Figure 2 .
Figure 2. Log-log plot between the aggregated variance of the SWS spatial series and the scale.A linear relationship indicated the presence of scale invariance and scaling laws for three selected dates.

Figure 4 .
Figure4.Mass exponents for selected soil water storage spatial series from surface to different soil layers (cumulative storage) at 20 cm increment down to 140 cm in 2008 for a range of q (−15 to 15 at 0.5 increments).The solid line is a linear reference created following the UM model ofSchertzer and Lovejoy (1987) passing through (q = 0).

Figure 5 .
Figure 5.The width of the multifractal spectrum (α max -α min value) for soil water storage at different depths (20 cm increment) for all measurements completed during the study period.

Figure 6 .
Figure 6.Multifractal spectra of soil water storage spatial series measured at each 20 cm soil layer down to 140 cm in 2008, 2010 and 2011 for a range of q (−15 to 15 at 0.5 increments).

Figure 7 .
Figure 7. Multifractal spectra of soil water storage spatial series from surface to different soil layers (cumulative storage) at 20 cm increment down to 140 cm in 2008, 2010, and 2011 for a range of q (−15 to 15 at 0.5 increments).

Figure 8 .
Figure 8.The information dimension (D 1 ) for soil water storage at different depths (20 cm increment) over the whole measurement period.

Figure 11 .
Figure 11.Multifractal spectra of joint distribution of SWS at 0-20 and 20-40 cm measured on 22 October 2008.Contour lines show the joint scaling dimensions of the SWS measurement series.

Figure 12 .
Figure 12.Conceptual schematics showing vegetation development over time, dominant water loss processes and the scaling behavior of soil water storage at different depths.The figure is developed based on field observations and scaling analysis.The scale of the figure is arbitrary.
cm for the20-40, 40-60,   60-80, 80-100, 100-120, and 120-140 cm layers, respectively.Average SWS for a single measurement varied from 3.40 to 7.16 cm.The highest average SWS for the surface layer was observed on 29 June 2011.The study area received a large amount of spring snowmelt (2010 received 642 mm, double the annual average precipitation) and rainfall during 2011 leading to the high SWS in the surface layer (Weather Canada historical report).The lowest average SWS for the surface layer was observed on 23 August 2008, which was one of the driest summers within the 5-year study period.The highest average SWS (on 29 June 2011) at the surface layer gradually decreased to 6.55 cm at the deepest layer and the lowest average SWS (on 23 August 2008) at the surface layer gradually increased to 5.28 cm at the 120-140 cm layer (Table

Table 1 .
Maximum, minimum, and average soil water storage (cm) at different depths (20 cm increment) over the whole measurement period.

Table 2 .
Correlation coefficients between joint multifractal indices (α and β) (n = 440) of the surface layer with those from subsurface layers at 20 cm intervals in 2008.