Joint multifractal analysis for three variables: Characterizing the effect of topography and soil texture on soil water storage

Geoderma 334 (2019) 15–23

Contents lists available at ScienceDirect

Geoderma journal homepage: www.elsevier.com/locate/geoderma

Joint multifractal analysis for three variables: Characterizing the effect of topography and soil texture on soil water storage

T

Asim Biswas School of Environmental Sciences, University of Guelph, Guelph, Ontario N1G 2W1, Canada

A R T I C LE I N FO

A B S T R A C T

Handling Editor: Morgan Cristine L.S.

Multifractal analysis describes heterogeneity in the distribution of a variable by characterizing and summarizing the variability across scales. Joint multifractal analysis has been widely employed to characterize scale relationships between two variables co-existing along a single geometric support. In this study, joint multifractal analysis was used for three variables coexisting in the same geometric support to describe the influence of topography (relative elevation) and soil texture (sand content) on water storage within a soil profile. Soil water storage, as well as sand content and relative elevation, were measured down to 1.4 m depth along a 576 m long transect in the hummocky landscape of central Saskatchewan, Canada. Joint multifractal analysis was conducted to consider both the strange attractor formalism and the method of moments. The variability in soil water storage, sand content and relative elevation was scale dependent and showed multifractal behavior. The spatial variability in relative elevation was strongly reflected on water storage across the analyzed spatial scales but the joint multifractal spectrum for sand content and water storage suggested a lower degree of correlation. The change in multifractality was also observed when there was high variability in relative elevation and texture. This clearly demonstrated the capability of joint multifractal analysis to completely characterize the scaling behavior among three variables.

Keywords: Power law Fractal theory Scaling Partition function Soil water

1. Introduction Information on spatial and temporal distribution of soil water is a key input in monitoring soil water balance and the global hydrological cycle, assessing land-atmospheric interactions, understanding a large number of surface and subsurface hydrological processes, testing the performance of various engineered covers, and validation of climate and hydrological models (Koster et al., 2004; Rodriguez-Iturbe et al., 1999; Sivapalan, 1992; Western et al., 2002). Many physical factors (e.g. topography, soil properties) and environmental processes (e.g. rainfall, evapotranspiration, runoff, and snow melt) operating at different intensities over a variety of scales (Entin et al., 2000; Goovaerts, 1998) give rise to complex and nested patterns or spatial (temporal) organization (Western et al., 1999) in soil water variability as a function of spatial (temporal) scales (Biswas and Si, 2011c; Gomez-Plaza et al., 2000; Kachanoski and Dejong, 1988). Basically, the spatial variability of soil water is the reflection of the spatial variability of hydrological processes. Thus, information on soil water variability at the scale of measurement provides indication on the underlying soil hydrological processes at that scale (Banerjee et al., 2011; Goovaerts, 1998). However, as the variability in soil water is controlled by several factors and processes operating at multiple scales, the multi-scale

nature of hydrological processes should be characterized and quantified to better understand the hydrological cycle (Biswas and Si, 2011b; Biswas and Si, 2011c; Gomez-Plaza et al., 2000). For example, at small catchment and hill slope scale, factors like water routing processes (Beven and Germann, 1982; Dunne and Black, 1970; Moore et al., 1988), differential radiation effects (Moore et al., 1993; Western et al., 1999), heterogeneity in soil (Famiglietti et al., 1998; Hu et al., 1997; Seyfried, 1998) and vegetation (Hupet and Vanclooster, 2002; Qiu et al., 2001) affect soil water content/storage patterns on and within the landscape. In contrast, atmospheric, geologic, and climatic variability determine the organization of soil water over large area (Brocca et al., 2007; Entin et al., 2000; Schneider et al., 2008; Seyfried, 1998). The contribution of certain processes, though dominant at one scale, may have weaker control at other scales and relatively smaller contribution towards overall hydrological dynamics. Now at a scale, the effect from weaker processes can be masked by other dominant processes (Biswas and Si, 2011a; Ji et al., 2016; Kachanoski and Dejong, 1988). Multi-scale variability of soil water and thus hydrological processes make hydrological studies challenging and managing decisions difficult at a scale other than the scale of measurement. Therefore, characterizing variations at the scale of measurement and transferring information from one scale to another, also known as scaling, is critical

E-mail address: [email protected]. https://doi.org/10.1016/j.geoderma.2018.07.035 Received 6 October 2017; Received in revised form 18 July 2018; Accepted 27 July 2018 0016-7061/ © 2018 Elsevier B.V. All rights reserved.

Geoderma 334 (2019) 15–23

A. Biswas

indices (Kravchenko et al., 2000; Zeleke and Si, 2004), grassland productivity and terrain indices (Banerjee et al., 2011), volume and number based soil particle size distribution (Li et al., 2011), and soil water retention parameters and soil texture (Wang et al., 2011). However, to date, joint multifractal analysis has been used to explore the relationship between only two variables at different scales. The scaling relationship between soil water storage and individual factors may not explore the relationship sufficiently as the overall effect from multiple factors could be highly nonlinear and non-additive (Biswas and Si, 2011b). For example, among other factors, soil water storage is known to be controlled by topography (Brocca et al., 2010; Western et al., 1999), soil texture (Cosh et al., 2008; Vachaud et al., 1985) or a combination (Biswas et al., 2012; Tallon and Si, 2004). In this situation, characterizing the joint variability between either soil water storage and topography or soil water storage and soil texture at multiple scales may not reveal the complete picture. Extending the joint multifractal analysis (Meneveau et al., 1990) to multiple variables (PavonDominguez et al., 2015) could provide a complete picture of soil water storage spatial variations as controlled by different factors at different scales. Therefore, the objective of this study was to extend the joint multifractal analysis for three variables to characterize the effect of topography (e.g. relative elevation) and soil texture (e.g. sand) on soil water storage at multiple scales. A 3D visualization of joint multifractal spectrum was used to describe the effect of multiple variables at multiple scales.

for understanding complex hydrological dynamics. Semivariogram, the central to geostatistical analysis and a secondorder statistical moment of a spatial variable (Journel and Huijbregts, 2003), has been used to characterize the spatial pattern in soil water (Brocca et al., 2007; Entin et al., 2000; Western et al., 2002). The underlying assumption of this is the ‘intrinsic hypothesis’ where all the regionalized variables are considered Gaussian and it only depends on the separation distance or lag distance between samples (stationary process) (Journel and Huijbregts, 2003). Nevertheless, it can effectively reveal the spatial distribution and autocorrelation of soil water within the study area. However, the complex behavior of soil water variability is created by both irregularity and structure for different length scales resulted in several closely-spaced and nested scales of variations, occurring over short distances (Burrough, 1983a; Goovaerts, 1997). To overcome the inadequacy of a single ‘range’ to account for abrupt changes of the mean of the target soil property e.g. soil water storage and the lack of deterministic similarity of variations at different scales, nested semivariogram model can be calculated (Burrough, 1983b). In this model, several non-overlapping and independent random functions each with its own weight (based on dominance and contribution) and scale were used to represent the underlying soil processes determining variability of soil water. However, practical application of this model requires a priori knowledge on the dominant scales and magnitudes of the processes controlling the pattern of soil water, which is rarely available (Zeleke and Si, 2006). Moreover, having quantified these relationships at multiple scales, it is required to transfer information from one scale to another for developing connection among ‘process scale’, ‘measurement scale’ and ‘management scale’ (Biswas and Si, 2011d) but limited by its capability. Additionally, traditional correlation analysis, commonly used to identify the dominant factors of soil water storage (Biswas et al., 2012; Vachaud et al., 1985), can quantify the relationship only at the scale of measurement and needs to be analyzed across scales. The scaling of a soil property e.g. soil water is possible if the distribution of some statistical parameters (e.g. variance) within a geographical support space or spatial domain remains the same at other geographical support spaces or spatial domains (Cheng, 1999). This means that the feature in the distribution of soil water will not change if the spatial domain is multiplied/divided by a common factor, also known as scale-invariance or self-similarity (Hu et al., 1997; Kim and Barros, 2002). Therefore, the probability of measuring a value will vary inversely as a power of that value and is known as power function, a typical scaling process. As the spatial distribution of soil water is known to follow the power law function (Hu et al., 1997; Ji et al., 2016; Kim and Barros, 2002; Mascaro et al., 2010), spatial variability can be investigated and characterized using fractal theory (Mandelbrot, 1982). As the soil water storage could be a response of complex nonlinear processes acting over a variety of scales, the spatial distribution characterization and quantification requires multiple scaling indices (multifractal scaling). Therefore, multifractal scaling indices can provide insights into the interrelationships between systems and the organization about the underlying mechanism (Cheng, 1999; Martin-Sotoca et al., 2018). For example, Ji et al. (2016) and Kim and Barros (2002) observed a multifractal behavior of surface soil water as a result of the temporal evolution of wetting and drying in a semi-arid and humid environment, respectively. Similarly, Mascaro et al. (2010) reported the multifractal behavior of soil water and used it to develop downscaling models for remote sensing measurements. As the scaling properties of the spatial patterns of soil water storage change with scales (Brocca et al., 2010; Ji et al., 2016; Kim and Barros, 2002), it is necessary to understand the relationship with controlling factors at different scales. If the multifractal behavior of soil water storage co-exist with other controlling factors in the same spatial domain, the relationship can be jointly analyzed using the joint multifractal analysis (Meneveau et al., 1990). Joint multifractal analysis has been used to characterize the interaction between crop yield and terrain

2. Materials and methods 2.1. Study site and data collection The data (soil water storage, relative elevation and sand) used in this study were collected along a north-south direction transect of 128 points separated by 4.5 m regular intervals. The transect was established in 2004 within the St. Denis National Wildlife Area (52°12′N lat. and 106°50′W long.) located in central Saskatchewan, Canada as part of a bigger project to study soil water dynamics in the Prairie Pothole Region. Several other publications have used part of the dataset collected over the years to answer various questions regarding soil water dynamics. More detailed information on the study site, actual data set, data collection procedures and publications using the dataset can be found in Biswas et al. (2012). Briefly, the landscape of the study area is mainly hummocky with a complex sequence of slopes (10 to 15%) extending from different sized rounded knolls to depressions. The depressions, also known as potholes, were formed from buried ice chunks during the last deglaciation (Huel, 2000) and almost act as an individual micro-watershed. This is a landscape characteristic of the North American Prairie pothole region covering approximately 780,000 km2 area from the north-central United States to south-central Canada (National Wetlands Working Group, 1997) and it is very important for its unique hydrological and ecological functions. The soil of the study area is mainly Dark Brown Chernozem developed from moderately fine to finely textured, moderately calcareous glacial till (Saskatchewan Centre for Soil Research, 1989). The climate is semi-arid with long-term annual average precipitation of 360 mm and an annual average air temperature of 2 °C (AES, 1997). The vegetation of the study area is mixed grass. The surface soil water (0–20 cm) was measured using a time domain reflectometry (TDR) probe and a metallic cable tester (Model 1502B, Tektronix, Beaverton, OR, USA). Sub-surface water content was measured at every 20-cm depth down to 140-cm using a neutron probe (Model CPN 501 DR Depthprobe, CPN International Inc., Martinez, CA, USA). A site-specific calibration was used for the neutron probe (Biswas et al., 2012) and a standard calibration equation (Topp and Reynolds, 1998) was used for the TDR measurements. Soil water storage (SWS) was calculated from soil water content and depth of each layer and added together to get the total SWS in the profile (0–140 cm). In this 16

Geoderma 334 (2019) 15–23

A. Biswas

properties can be determined. For example, effect of some nonlinear instantaneous variations in soil properties on the variability at different scales can be quantified following multifractal theory. Therefore, spatial random variable can have multifractal properties for any q with ∞ ≤ q ≤ −∞. Similarly, the curvature and linearity/nonlinearity of the τ(q) is often used to measure the degree of multifractality. Characterizing multifractal behavior can provide new insight on the spatial pattern in soil properties (e.g. soil water). Traditional and ordinary statistics including semivariogram, lacunarity analysis, autocorrelation is associated with local multifractality around the mean and help evaluate the local properties around the mean. In soil science application, where extreme values are often of interest, the entire multifractal spectrum and statistics based on higher order statistical moments provide improved information on the variability of soil properties and the basis for transferring information from one scale to another. One such strong application is downscaling information from satellite data (Kim and Barros, 2002; Mascaro et al., 2010). While multifractal analysis has been used to characterize the scaling property of one variable, joint multifractal analysis has been used to characterize the joint distribution of scaling properties of two variables coexisting in the same geometric support. This study extends the joint multifractal analysis to three variables using both the strange attractor formalism (Grassberger, 1983; Halsey et al., 1987; Hentschel and Procaccia, 1983) and the method of moments (Evertsz and Mandelbrot, 1992). The overall procedure was subdivided into five steps and explained using the examples of three spatial series of soil water storage (SWS), relative elevation (RE) and sand (Sand). Step 1- Subdivision of the spatial series: The spatial domain of the data series was divided into a smaller number of segments based on a rule that generated self-similar segments. In this study, the transect was divided into nini non-overlapping segment of the spatial resolution equal to the sampling interval (δini = 20 = 4.5 m = 1 data). Therefore, the spatial resolution of the initial segment, j for any series [SWS]ini, [RE]ini or [Sand]ini was equal to the sampling interval of 4.5 m. The transect was subsequently subdivided into n non-overlapping segments for a spatial resolution ranging from δini to δfinal = 27 = 128 sampling points (×4.5 m sampling interval) = 576 m. For each spatial resolution (δ), the sum of the values of [SWS]ini, [RE]ini and [Sand]ini within that segment was assigned to the interval. Step 2- Calculation of probability mass function: For the multifractal measure, the probability mass function for each segment, ci, was defined as (Kravchenko et al., 2000; Zeleke and Si, 2004):

study, SWS was measured on 17 July 2007 and used as an example. The topographic survey of the study site was completed using Light Detection and Ranging (LiDAR) survey at 5 m resolution and geographic location of the sample points and relative elevation (RE) was recorded using a Trimble Pro XRS Global Positioning System (Trimble Navigation, Sunnyvale, CA, USA). Surface soil texture was determined using a HORIBA LA-950 Laser diffraction particle-size analyzer (HORIBA Scientific Inc., Edison, NJ, USA). 2.2. Data analysis To describe complexity and self-similarity in nature, various fields of science have adopted and used the concept of fractals and multifractals, which are a consequence of self-similarity arising from scaleindependent processes. The fractal concept was introduced to represent irregular geometry by its Hausdorff dimension, which is greater than its topological dimension (Mandelbrot, 1982), while the multifractal describes the measures from which spatially intertwined fractals, each characterized by singularity strength and fractal dimension, can be derived (Cheng, 1999; Chhabra et al., 1989; Evertsz and Mandelbrot, 1992; Halsey et al., 1987; Meneveau et al., 1990; Schertzer and Lovejoy, 1991). It is well known that the fractal property is determined by power-law relationship between the measure (e.g. soil water storage) and the measuring unit (e.g. spatial distance or sampling intervals). The self-similarity or self-affinity of the fractals as characterized by the scale-invariant property of the power-law function (changing the measuring unit does not change the function type) describes multifractal behavior of the property of interest (Cheng, 1999; Cheng and Agterberg, 1996; Feder, 1988; Schertzer and Lovejoy, 1991). Various complex processes such as multiplicative cascade processes, diffusionlimited aggregation, turbulence and Brownian motions have been suggested as multifractal measures (Cheng, 1999; Evertsz and Mandelbrot, 1992; Feder, 1988; Schertzer and Lovejoy, 1991). Spatially intertwined multifractals (continuous or discrete) provides information on spatially related fractals or self-organized fractals and can describe the spatial variability of the measure (e.g. soil water) as a function of scale (Cheng, 1999; Cheng and Agterberg, 1996). Multifractal analysis has been used to study the spatial distribution of physical and chemical quantities with geometric support consisting of a set of points of spatial objects (Cheng and Agterberg, 1996). Various methods have been used in characterizing multifractal behavior and the details are beyond the scope of this paper. Most often three functions are used; τ(q), called ‘mass exponent’ function, α(q), Coarse Holder exponent and f(α), fractal/multifractal spectrum and can be related to explain the spatial variability. For example, the values of τ(q) at q = 0, 1, and 2, and the local properties of τ(q) at q = 1 is useful in describing and characterizing conventional second-order moment statistics (Cheng, 1999; Schertzer and Lovejoy, 1991). The second-order multifractality for q = 1 and 2 is often sufficient to derive exact and appropriate expressions for the covariance, autocorrelation and semivariogram function of continuous random variables in spatial statistics and thus can provide similar information. For example, semivariogram function is primality determined by the value τ(2). Similarity lacunarity, a measure the deviation of a geometric object such as fractal from translational invariance, can be quantified using τ(q) for q ranging 0 and 2. For example, geometric objects with relatively low lacunarity are homogeneous and translationally invariant. These statistical parameters are determined by the exponents of the first and second order statistical moments and known as local multifractality of the mean measure (Cheng, 1999). For some properties with highly skewed distribution or properties are dominated by relatively few high or low values, often these statistical moments are not enough to express the system. Transformation of data can be a choice but often restricted to the incomplete expression of variability. Now as the q values tend to be ∞ or −∞, the effect of the largest or the smallest values, respectively on the variability quantification of a highly skewed distributed

ci [SWS ](δ ) =

ci [RE ](δ ) =

[SWS ]i n ∑ j =ini1 ([SWS ]ini )j

[RE ]i n ∑ j =ini1 ([RE ]ini )j

ci [Sand](δ ) =

[Sand]i n ∑ j =ini1 ([Sand]ini )j

(1)

(2)

(3)

where [SWS]i, [RE]i and [Sand]i were the sum of values of soil water storage, relative elevation and sand within the interval i for a specific n n spatial resolution of δ and ∑ j =ini1 ([SWS ]ini )j , ∑ j =ini1 ([RE ]ini )j and n

∑ j =ini1 ([Sand]ini )j were the sum of all values of the spatial series of soil water storage, relative elevation and sand, respectively. Therefore, the probability mass function (ci[SWS](δ), ci[RE](δ), ci[Sand](δ)) was the ratio of the relative weight or density of each segment to the sum of all values of the spatial series. Step 3- Calculation of joint partition function: Following the methods of moments, the joint partition function of two variables e.g. SWS and RE was calculated as (Evertsz and Mandelbrot, 1992): χ (q[SWS] , q[RE ] , δ ) = 17

[ci [SWS] (δ )]q[SWS] [ci [RE ] (δ )]q[RE ] n ∑i = 1

{[ci [SWS] (δ )]q[SWS] [ci [RE ] (δ )]q[RE ] }

(4)

Geoderma 334 (2019) 15–23

A. Biswas

where n indicated the number of segments for a specific δ value. The exponents (real numbers; q[SWS], q[RE]) were then used to weigh the joint distribution of SWS and RE. The current work extended and corroborated the existence of scaling properties in the joint distribution of three variables. The joint partition function of three variables (SWS, RE and Sand) depended on the values of the three exponents and was calculated for each segment, δ as follows:

The fractal dimension, f(α[SWS], α[RE], α[Sand]) can be defined from the scaling relationship

N (q[SWS] , q[RE ] , q[Sand] , δ ) ≈ δ−f (α[SWS] , α[RE ] , α[Sand] )

where N(q[SWS], q[RE], q[Sand], δ) was the number of intervals of a given spatial resolution δ. For a given combination of α values, the f (α[SWS], α[RE], α[Sand]) was the set of intervals that corresponded to a combination of the singularities α[SWS], α[RE] and α[Sand] and was calculated from (Chhabra and Jensen, 1989; Chhabra et al., 1989; Meneveau et al., 1990):

χ (q[SWS] , q[RE ] , q[Sand] , δ ) =

[ci [SWS] (δ )]q[SWS] [ci [RE ] (δ )]q[RE ] [ci [Sand] (δ )]q[Sand] n ∑i = 1

{[ci [SWS] (δ )]q[SWS] [ci [RE ] (δ )]q[RE ] [ci [Sand] (δ )]q[Sand] }

(5)

f (α[SWS] , α[RE ] , α[Sand] ) = q[SWS] α[SWS] + q[RE ] α[RE ] + q[Sand] α[Sand]

The value of three exponents (q[SWS], q[RE], q[Sand]) ranging from −5 to 5 at 0.25 increments were selected for this study. Positive value of qi magnifies the effect of large numbers and diminishes the effect of small numbers while for negative values of qi magnifies the effect of small numbers and diminishes the effect of large numbers in the spatial series i. The q values in a similar term can be explained as the order of statistical moments and be used to characterize the variability in a soil property by examining the effect of the magnitude of the data. Similarly, the effect of one or more variables on the others can also be characterized by varying the q values. For a special case of q = 0 for two variables, the joint partition function is simplified and resembles the partition function of the third variable and can be used to characterize the spatial variability. For the monofractal and multifractal measures, the scaling property of the joint partition function was written as:

χ (q[SWS] , q[RE ] , q[Sand] , δ ) ≈ δ τ (q[SWS] , q[RE ] , q[Sand] )

− τ (q[SWS] , q[RE ] , q[Sand] )

α[RE ] (q[SWS] , q[RE ] , q[Sand] ) =

(11)

The multifractal spectrum for a spatial series reaches a maximum value of f(α[SWS], α[RE], α[Sand]) equal to 1, which corresponded to the capacity dimension of the geometric support for the spatial series. Now this joint multifractal spectrum for three variables could be transformed to a joint multifractal spectrum for two variables when the statistical exponent of one variable is set to zero and multifractal spectrum for one variable when the statistical exponents for two variables are set to zero. A general form of multifractal spectrum can be presented as f(αi) versus αi(q) for spatial series i and could be used to study the effect of one or more variables on the spatial distribution of the rest of the variables. In this study, I examined the effect of such environmental factors as topography (e.g. relative elevation) and soil texture (e.g. sand) on soil water storage. The analysis was conducted in MATLAB (The MathWorks Inc., Natrick, MA, USA. Release: 2013a.) after modifying the pseudo-code programming language provided by PavonDominguez et al. (2015).

(6)

where τ(q[SWS], q[RE], q[Sand]) was the mass exponent function and depended on the exponents q[SWS], q[RE] and q[Sand]. Step 4- Calculation of joint mass exponent function: The mass exponent is related to the probability of the mass function in each segment to the size of the segment or spatial resolution. For a single variable, a linear relationship between the log-log plot of the partition function and the spatial resolution, δ, also known as mass exponent function, can be established when the analyzed measures show either monofractal, or multifractal, behavior. Similarly, in this study, the joint mass exponent function was estimated from the slope of the linear fits of the log-log plot of the joint partition function and the spatial resolution and was dependent on the three exponents. This dependence of the mass exponent function against the exponents confirms the existence of the multifractal nature, while linear fits from the logarithmic plots of the joint partition function and the spatial resolution reveal the range of spatial scales δmin – δmax, for which the multifractal nature was exhibited. Step 5- Calculation of joint multifractal spectrum: The joint multifractal spectrum characterized the joint distribution of the local exponents α[SWS], α[RE] and α[Sand] and was presented as f (α[SWS], α[RE], α[Sand]) versus α(q[SWS], q[RE], q[Sand]) (Meneveau et al., 1990). The Lipschitz-Hölder or singularity exponents α(q[SWS], q[RE], q[Sand]) quantified the strength of the measures of singularities or local scaling index. The singularity exponents were determined after the Legendre transformation (Evertsz and Mandelbrot, 1992) of the mass exponents for each statistical moment and were expressed as:

α[SWS] (q[SWS] , q[RE ] , q[Sand] ) =

(10)

3. Results and discussion 3.1. Spatial pattern of soil properties An average of 40.98 cm with a minimum of 23.95 cm and maximum of 68.08 cm soil water was stored in the profile along the transect (Fig. 1). The broad range (44.13 cm) was due to the variable water storage in the soil profile at different landscape positions. Late snowmelt (mid-May) due to extended winter in the study area along with heavier early summer precipitation created a large variation in SWS (Biswas et al., 2012). Though the total precipitation in 2007 (366 mm) was very similar to the average of long-term precipitation (360 mm) in the area, concentration of rainfall events in early summer contributed a

dτ (q[SWS] , q[RE ] , q[Sand] ) dq[SWS]

(7)

dτ (q[SWS] , q[RE ] , q[Sand] )

α[Sand] (q[SWS] , q[RE ] , q[Sand] ) =

dq[RE ]

(8) Fig. 1. Spatial distribution of soil water storage (top), relative elevation (middle) and sand content (bottom) along the transect. Histograms at the right side of the spatial distribution figure are showing the statistical distribution of the measurements along with theoretical normal distribution.

dτ (q[SWS] , q[RE ] , q[Sand] ) dq[Sand]

(9) 18

Geoderma 334 (2019) 15–23

A. Biswas

low correlation (Biswas and Si, 2011a; Biswas and Si, 2011b). The negative association at the scale of measurement between the Sand and SWS can be explained as the low water holding capacity of sandy soil (Fig. 2) (Famiglietti et al., 1998). However, this relationship might also be different at different scales.

large amount of water to the landscape. The depression along the transect was completely saturated and had standing water during measurement. For example, two depressions located around 120 m and 230 m from the origin of the transect (Fig. 1) had standing water and the SWS at those locations were calculated based on the bulk density assuming complete saturation. However, the knolls and the long slope at the end of the transect stored a comparatively small amount of water, creating a broad range in the measurements (Biswas et al., 2012; Biswas and Si, 2011a). This also created a large standard deviation (SD) (0.78 cm) in the measurements with the coefficient of variation (CV) of 22%. However, the measurements exhibited near-normal distribution (Fig. 1) with a skewness and kurtosis value of 0.69 and 0.25, respectively. The average RE along the transect was 3.99 m with a minimum of 2.39 m and a maximum of 5.31 m (Fig. 1). The SD of the measurement was 0.76 m and CV was 19%. The average sand content along the transect was about 30% with a range between 11 and 55% (Fig. 1). A little larger SD (11%) of measurement was observed for sand content with the CV of 37%. Both RE and Sand exhibited a near-normal distribution. The knolls stored a small amount of water compared to the depressions and created a spatial pattern in SWS that was almost like a mirror image of the spatial pattern of RE (Biswas and Si, 2011b). However, the correlation coefficient (r2 = 0.06) was not able to identify the association between them (Fig. 2). The Pearson correlation analysis explained the global association between SWS and RE irrespective of scales. The relationship at the scale of measurement might have been different from other scales and the positive and negative relationship at different scales might have neutralized each other exhibiting an overall

3.2. Joint multifractal analysis The integrated scaling relationship among SWS, RE and Sand was examined over a range of statistical moments or q exponents by exploring the scaling dependence between their joint partition function and the spatial resolution. Linear fits were obtained (linear regression coefficient of determination, R2 > 0.999) in the log-log scatter plot between the joint partition function and spatial resolution between δini = 4.5 m (single sampling interval) to δfinal = 576 m (length of transect). The linear relationship indicated the presence of (multi-) fractal behavior of the joint distribution of SWS, RE and Sand (Brocca et al., 2010; Ji et al., 2016; Mascaro et al., 2010). In reality, there will be one fractal dimension (one point in the spectrum) for a spatial series that exhibits a perfect monofractal distribution. This means, that the behavior of SWS measured at different spatial resolutions will only be a function of one scaling exponent and the information will be a simple multiplication between the information at the smallest spatial resolution and the scaling exponent. The joint multifractal spectrum of multiple variables each with perfect monofractal distribution will also have a single point. However, if the fractal dimension changes with the change of moments or q values, there will be multiple points in the spectrum and the degree of heterogeneity or spatial variability can be studied from the range of fractal dimension values. This means, as we change the moment of statistical parameters, the spatial relationship will not be a simple multiplication with information at the smallest resolution as the exponent will change. As the exponent changes, it creates a nonlinear pattern in the multifractal spectrum, f(α) and the degree of heterogeneity increases as the range between αmax and αmin increases, and so does the spatial variability. The spatial variability of SWS in relation to RE and Sand was quantified in this study. The joint multifractal spectrum of these variables can be presented as a plot of fractal dimension f(α[SWS], α[RE], α[Sand]) for the set of segments corresponding to each combination of singularity exponents; α[SWS], α[RE] and α[Sand] (Fig. 3a). The spectrum is the set of surfaces at different q values ranging from −5 to 5 at intervals of 0.25. For a q value of any variable (e.g. SWS), one can separate a slice of the surface from the multifractal spectrum that contains specific characteristics for that variable (e.g. SWS) in relation to two other variables (e.g. RE and Sand). The relationship between SWS and the RE and Sand was established based on the shape of the surface. The joint multifractal spectrum reached the maximum value for f(α[SWS], α[RE], α[Sand]) equal to 1 when the statistical exponents for all three variables (SWS, RE and Sand) were equal to zero (q[SWS], q[RE], q[Sand] = 0). This maximum value is known as the capacity dimension of the spatial support (Zeleke and Si, 2006). First, the joint multifractal relationship between SWS and RE was analyzed using certain fixed effects of Sand. For example, three specific q[Sand] exponents were selected to study the effect of RE on SWS distribution; q[Sand] = 5 and q[Sand] = −5 representing the situation with high and low Sand content, respectively and q[Sand] = 0 representing average Sand content. For each q[Sand] value, q[SWS] and q[RE] created a 2D surface. Fig. 3b depicted the three surfaces generated for q[Sand] = −5, 0 and 5 in relation to the overall joint multifractal spectrum while Fig. 3c depicted only three surfaces. These surfaces were further projected on the q[SWS] - q[RE] plane to study the effect of RE on SWS at fixed Sand effects (Figs. 4a–c). For low Sand content (q[Sand] = −5), the joint distribution of SWS and RE showed a widening spectrum from top left to the bottom right (Figs. 4a and 5). On the bottom right part of the spectra, high values of α[SWS] and low values of α[RE] were found. This indicated that the high SWS were related to the

Fig. 2. Scatter plot of soil water storage and relative elevation (bottom) and soil water storage and sand (top). The black solid line indicated the regression line with correlation coefficient (r) and the probability (p) of statistical significance given within each figure. 19

Geoderma 334 (2019) 15–23

A. Biswas

Fig. 3. a) Joint multifractal spectrum of the spatial series of soil water storage (SWS), relative elevation (RE) and sand (Sand). The surfaces corresponding to the q[RE] = −5, 0 and 5 projected on the joint multifractal spectrum (b) and separated individually (c) and the surfaces corresponding to the q[Sand] = −5, 0 and 5 projected on the joint multifractal spectrum (d) and separated individually (e).

the depressions, longer infiltration time makes the SWS less variable. Similarly, at the higher sand content (q[Sand] = 5), the effect of low (q[RE] = −5) and high (q[RE] = 5) RE was also analyzed and presented as solid black symbols in Fig. 4c. For both cases, smaller values of f (α[SWS], α[RE], α[Sand]) indicated a weaker variability in SWS in either depressions or on knolls when there is higher sand content (Biswas et al., 2012; Biswas and Si, 2011b; Western et al., 1999). This weaker effect may be due to faster gravity-driven infiltration in sandy soil irrespective of the topographic positions making SWS relatively less variable across scales (Famiglietti et al., 1998; Ji et al., 2016; Seyfried, 1998). This clearly indicated a stronger association between SWS, elevation and soil texture (Biswas et al., 2012) at different scales. The effect of sand content on SWS at fixed RE was also studied (Figs. 4d–f). The variability in SWS was also prominent at varying sand content at different landscape positions or at fixed RE. This was visible from the larger differences of f(α[SWS], α[RE], α[Sand]) values for different sand content at a fixed RE. For example, at depressions or at low RE (q[RE] = −5), the effect of small amounts of sand (q[SWS] = −5) was quite different than the effect of large amounts of sand (q[SWS] = 5) (Figs. 4d and 5). The large difference in the maximum and minimum values of α[SWS] indicated the scaling heterogeneity in SWS resulting from smaller amounts of sand or larger amounts of silt and clay. However, the effect was less prominent when there was a large amount of sand. Large amounts of sand in the depressions affected SWS from the higher infiltration time and faster infiltration rate (Biswas et al., 2012; Famiglietti et al., 1998). On the other hand, higher amounts of clay and

low RE at low sand content (Biswas et al., 2012). In hummocky landscapes, water redistributes from knoll to depressions. Similarly, water erosion is an issue where smaller sized particles get redistributed from knolls to depressions leading to high SWS in depressions with low sand or high amount of clay (Biswas et al., 2012). On the top left part of the spectra, a relationship between low SWS and high RE was also established (low α[SWS] and high α[RE]). Therefore, the well-known direct statistical relationship between SWS and RE was maintained across all analyzed spatial scales (Ji et al., 2016). Moreover, the multifractal spectra, f(α[SWS], α[RE], α[Sand]) value moves from the highest to the lowest indicating a greater heterogeneity in the space interval covering high RE and low SWS. For low Sand content (q[Sand] = −5), the effect of low (q[RE] = −5) and high (q[RE] = 5) RE was also analyzed and presented as solid black symbols in Fig. 4a. For low Sand content (q[Sand] = −5) and low RE (q[RE] = −5), the values of f (α[SWS], α[RE], α[Sand]) slowly decreased from low α[SWS] to high α[SWS] while for low Sand and high RE (q[RE] = 5), the values of f (α[SWS], α[RE], α[Sand]) showed a comparatively faster decrease (Figs. 4a and 5). At the same time, the large range of f(α[SWS], α[RE], α[Sand]) values at low Sand and high RE indicated greater heterogeneity compared to low Sand and low RE. This means that the depressions with small amounts of Sand (or high amounts of silt and clay) showed weaker variability is SWS compared to the knolls with low Sand content. This persisted over scales. This may be due to the reduced infiltration time on knolls with higher silt and clay content leading to surface runoff to depressions (Famiglietti et al., 1998; Groenendyk et al., 2015). While in 20

Geoderma 334 (2019) 15–23

A. Biswas

Fig. 4. Joint multifractal spectra between soil water storage (SWS) and relative elevation (RE) for q[Sand] = −5 (a), 0 (b) and 5 (c) and the joint multifractal spectra between SWS and sand (Sand) for q[RE] = −5 (d), 0 (e) and 5 (f).

across scales at low sand content on knolls. This means that simple scaling indices can explain the variability of SWS as affected by sand on knolls (Brocca et al., 2010; Ji et al., 2016; Kim and Barros, 2002; Mascaro et al., 2010). A specific case of the joint multifractal spectra is the multifractal spectrum of SWS that can be obtained geometrically from the intersection of the surface for q[Sand] = 0 and q[RE] = 0 (Fig. 4b and c). This corresponded to the scaling behavior of SWS as a single variable (Fig. 5) without the effect of soil texture and topography. Additionally, the influence of other variables can also be avoided by calculating a single spectrum from the multifractal analysis of the SWS spatial series as a single variable (Fig. 6). As the whole spatial series of SWS was considered here, f(α)max was equal to 1. However, as other factors such as soil texture and topography came into action, the maximum value did not reach unity due to interactive influences. For example, at high or low sand content at depressions or on knolls, the maximum value of f (α[SWS], α[RE], α[Sand]) did not reach unity and showed a strongly marked non-symmetric shape with right tails sometimes being longer than left ones (knolls with high or low sand content) or left tails being longer than rights ones (depression with high or low sand content) (Fig. 5). For longer right tails, the α[SWS] reached 1.026 at f(α) of 0.58 and was related to the low SWS at the knolls with small amounts of sand content or large amounts of clay and silt. However, the α[SWS] of 1.073 at f(α) of 0.60 was related to the knolls with large amounts of sand. Generally, the high α[SWS] value indicated the presence of low SWS, the low f(α) values suggested that the scales covering low SWS was rather small and localized and the higher variability was observed for the low SWS spatial distribution (Brocca et al., 2010; Ji et al., 2016; Mascaro et al.,

Fig. 5. Joint multifractal spectra (f(α) vs α) between the spatial series of soil water storage (SWS) and relative elevation (RE) and sand (Sand) for a combination of selected moment orders (q).

silt will have a different influence on infiltration and thus, SWS (Figs. 4d and 5). Though similar values of α[Sand] were observed at knolls (or higher topographic positions with high RE values), the difference in variability at small and large amounts of sand was not very prominent (Figs. 4f and 5). This indicated that sand had a less prominent effect on SWS on knolls. However, unlike all these situations, the SWS spatial series showed an almost monofractal type of behavior 21

Geoderma 334 (2019) 15–23

A. Biswas

and SWS and RE were observed with varying degrees of multifractality. This means that the spatial variability of RE and Sand was strongly reflected on the SWS across the studied spatial scales. The effect of RE seemed to be relatively stronger on SWS than Sand. Depressions with smaller amounts of sand seemed to affect SWS strongly and vary across scales over depressions with a large amount of sand. On the contrary, knolls with large amounts of sand seemed to have a stronger influence on SWS scaling than with small amounts of sand which may be associated with infiltration rate and time. Therefore, unlike correlation analysis, the joint multifractal spectra better reflected the variability of SWS as affected the RE and sand content. It supplements and completes the information provided by statistics as the relationships were examined over a wider range of spatial scales and exponent orders or moments. This provides an important information on soil spatial variability for data with extreme distributions. Finally, this study provided the first step to analyze and quantify the joint scaling properties of more than two variables. Future work with more examples from different fields and with more than three variables would make the methodology robust and help explain the variability and underlying processes that are controlled by multiple variables coexisting in the same geometric support.

Fig. 6. Multifractal spectra (f(α) vs α) of spatial series of soil water storage (SWS), relative elevation (RE) and sand (Sand) for a range of moment orders (q) from −5 to 5 with an increment of 0.25 moment.

2010). This was stronger at knolls with large amounts of sand and may be due to the effect of varied solar radiation (Biswas et al., 2012). On the contrary, the scale covering high SWS was significant as indicated by the high f(α) values at the left side of the spectrum (Fig. 5). An opposite pattern in the f(α) values was observed in depressions with larger or smaller amounts of sand. The longer left tails with smaller α[SWS] and higher f(α) values were related to the high values of SWS. Low α[SWS] indicated the presence of little SWS and higher f(α) suggested the scales covering large SWS was rather large (Fig. 5). However, in none of the cases did the maximum value of f(α) reach unity because the analysis was focused on the response of certain values of SWS (based on the q values selected) and not the entire spatial series. This suggested an interactive effect from soil texture and topography which was otherwise absent in single multifractal spectrum (Biswas et al., 2012; Biswas and Si, 2011b; Brocca et al., 2010) (Fig. 6). While the single multifractal spectrum mostly showed a symmetric structure with small variations in tail length mainly indicating a multifractal nature with regular distribution of scaling indices. However, the highest level of multifractality was observed for sand due to inherent variability of soil texture in glacial till materials of the hummocky landscape (Biswas et al., 2012; Ji et al., 2016). The RE and SWS showed an almost mirror image of the multifractal spectrum for each other; this may be due to the resemblance of the spatial patterns with depressions storing more water compared to knolls (Biswas et al., 2012). The opposite behavior of the tail length also supported the mirror image behavior.

Acknowledgements The dataset used in this study was collected under the supervision of Dr. Bingcheng Si and with the help of Dr. Henry Chau. I am grateful for their help and support. The dataset can be available on request. I am also grateful to Dr. Ana-Maria Tarquis for her thoughtful comments and discussion on the topic and the thought-provoking and intriguing comments from the reviewer and the Associate Editor (Christine Morgan). The funding for this work was provided by NSERC Discovery grant (RGPIN-2014-04100). I am also grateful to Dr. Pavon-Dominguez from Spain for the code and is available at https://link.springer.com/ article/10.1007/s00477-014-0973-5. References AES, 1997. Canadian Daily Climate Data for Western Canada. Environment Canada, Downsview, ON. Banerjee, S., He, Y., Guo, X., Si, B.C., 2011. Spatial relationships between leaf area index and topographic factors in a semiarid grassland: joint multifractal analysis. Aust. J. Crop. Sci. 5 (6), 756–763. Beven, K., Germann, P., 1982. Macropores and water flow in soils. Water Resour. Res. 18 (5), 1311–1325. Biswas, A., Si, B.C., 2011a. Identifying scale specific controls of soil water storage in a hummocky landscape using wavelet coherency. Geoderma 165 (1), 50–59. Biswas, A., Si, B.C., 2011b. Revealing the Controls of Soil Water Storage at Different Scales in a Hummocky Landscape. Soil Sci. Soc. Am. J. 75 (4), 1295–1306. Biswas, A., Si, B.C., 2011c. Scales and locations of time stability of soil water storage in a hummocky landscape. J. Hydrol. 408 (1–2), 100–112. Biswas, A., Si, B.C., 2011d. Scaling soil physical properties. In: Gliński, J., Horabik, J., Lipiec, J. (Eds.), Encyclopaedia of Agrophysics. Springer, The Netherlands, pp. 725–729. Biswas, A., Chau, H.W., Bedard-Haughn, A.K., Si, B.C., 2012. Factors controlling soil water storage in the hummocky landscape of the Prairie Pothole Region of North America. Can. J. Soil Sci. 92 (4), 649–663. Brocca, L., Morbidelli, R., Melone, F., Moramarco, T., 2007. Soil moisture spatial variability in experimental areas of central Italy. J. Hydrol. 333 (2–4), 356–373. Brocca, L., Melone, F., Moramarco, T., Morbidelli, R., 2010. Spatial-temporal variability of soil moisture and its estimation across scales. Water Resour. Res. 46. Burrough, P.A., 1983a. Multiscale sources of spatial variation in soil. I. The application of fractal concepts to nested levels of soil variation. J. Soil Sci. 34 (3), 577–597. Burrough, P.A., 1983b. Multiscale sources of spatial variation in soil. II. A non-Brownian fractal model and its application in soil survey. J. Soil Sci. 34 (3), 599–620. Cheng, Q.M., 1999. Multifractality and spatial statistics. Comput. Geosci. 25 (9), 949–961. Cheng, Q.M., Agterberg, F.P., 1996. Multifractal modeling and spatial statistics. Math. Geol. 28 (1), 1–16. Chhabra, A.B., Jensen, R.V., 1989. Direct determination of the f(α) singularity spectrum. Phys. Rev. Lett. 62 (12), 1327–1330. Chhabra, A.B., Meneveau, C., Jensen, R.V., Sreenivasan, K.R., 1989. Direct determination of the f (alpha) singularity spectrum and its application to fully-developed turbulence. Phys. Rev. A 40 (9), 5284–5294. Cosh, M.H., Jackson, T.J., Moran, S., Bindlish, R., 2008. Temporal persistence and

4. Summary and conclusions While traditional correlation analysis is the most common method to identify the dominant factor of SWS, the relationship depends on the size or length of the measurement scale and can explain the association with numerous factors only at the scale measurement. Multifractal scaling behavior of a SWS spatial series can be characterized based on the determination of scaling properties across scales through a robustly synthesized multifractal spectrum. The traditional joint multifractal spectrum can explain scaling behavior between two variables while multiple variables can simultaneously affect the distribution of SWS. This study extended the traditional joint multifractal analysis to examine the joint scaling behavior among three variables. An example SWS spatial series and the spatial series of two dominant controlling factors (RE and Sand) were considered in this study to implement the method and describe their joint scaling relationship. The scale dependence and multifractal nature of the studied SWS, RE, and Sand were confirmed and their joint influence was characterized over a range of scales. A strong and direct scaling relationship between SWS and Sand 22

Geoderma 334 (2019) 15–23

A. Biswas

Martin-Sotoca, J.J., Saa-Requejo, A., Borondo, J., Tarquis, A.M., 2018. Singularity maps applied to a vegetation index. Biosyst. Eng. 168, 42–53. Mascaro, G., Vivoni, E.R., Deidda, R., 2010. Downscaling soil moisture in the southern Great Plains through a calibrated multifractal model for land surface modeling applications. Water Resour. Res. 46 (8), W08546. Meneveau, C., Sreenivasan, K.R., Kailasnath, P., Fan, M.S., 1990. Joint multifractal measures: Theory and applications to turbulence. Phys. Rev. A 41 (2), 894–913. Moore, I.D., Oloughlin, E.M., Burch, G.J., 1988. A contour based topographic model for hydrological and ecological applications. Earth Surf. Process. Landf. 13 (4), 305–320. Moore, I.D., Gessler, P.E., Nielsen, G.A., Peterson, G.A., 1993. Soil attribute prediction using terrain analysis. Soil Sci. Soc. Am. J. 57 (6), 443–452. National Wetlands Working Group, 1997. The Canadian Wetland Classification System. University of Waterloo, ON. Pavon-Dominguez, P., Jimenez-Hornero, F.J., Gutierrez de Rave, E., 2015. Joint multifractal analysis of the influence of temperature and nitrogen dioxide on tropospheric ozone. Stoch. Env. Res. Risk A. 29 (7), 1881–1889. Qiu, Y., Fu, B.J., Wang, J., Chen, L.D., 2001. Soil moisture variation in relation to topography and land use in a hillslope catchment of the Loess Plateau, China. J. Hydrol. 240 (3–4), 243–263. Rodriguez-Iturbe, I., D'Odorico, P., Porporato, A., Ridolfi, L., 1999. On the spatial and temporal links between vegetation, climate, and soil moisture. Water Resour. Res. 35 (12), 3709–3722. Saskatchewan Centre For Soil Research, 1989. Rural Municipality of Grant, Number 372: Preliminary Soil Map and Report. Publication# SK372. Saskatchewan Centre Soil Research, University of Saskatchewan, Saskatoon, SK. Schertzer, D., Lovejoy, S., 1991. Non-Linear Variability in Geophysics: Scaling and Fractals. Springer Netherlands, The Netherlands. Schneider, K., Huisman, J.A., Breuer, L., Zhao, Y., Frede, H.G., 2008. Temporal stability of soil moisture in various semi-arid steppe ecosystems and its application in remote sensing. J. Hydrol. 359 (1–2), 16–29. Seyfried, M., 1998. Spatial variability constraints to modeling soil water at different scales. Geoderma 85 (2–3), 231–254. Sivapalan, M., 1992. Scaling of hydrologic parameterizations, 1. In: Simple Models for the Scaling of Hydrologic State Variables, Examples and a Case Study. Center for Water Research, University of Western Australia, Nedlands, WA, Australia. Tallon, L.K., Si, B.C., 2004. Representative soil water benchmarking for environmental monitoring. J. Environ. Inf. 4 (1), 28–36. Topp, G.C., Reynolds, W.D., 1998. Time domain reflectometry: a seminal technique for measuring mass and energy in soil. Soil Tillage Res. 47 (1–2), 125–132. Vachaud, G., Desilans, A.P., Balabanis, P., Vauclin, M., 1985. Temporal stability of spatially measured soil-water probability density-function. Soil Sci. Soc. Am. J. 49 (4), 822–828. Wang, Z.-Y., Shu, Q.-S., Xie, L.-Y., Liu, Z.-X., Si, B.C., 2011. Joint multifractal analysis of scaling relationships between soil water-retention parameters and soil texture. Pedosphere 21 (3), 373–379. Western, A.W., Grayson, R.B., Bloschl, G., Willgoose, G.R., McMahon, T.A., 1999. Observed spatial organization of soil moisture and its relation to terrain indices. Water Resour. Res. 35 (3), 797–810. Western, A.W., Grayson, R.B., Bloschl, G., 2002. Scaling of soil moisture: a hydrologic perspective. Annu. Rev. Earth Planet. Sci. 30, 149–180. Zeleke, T.B., Si, B.C., 2004. Scaling properties of topographic indices and crop yield: multifractal and joint multifractal approaches. Agron. J. 96 (4), 1082–1090. Zeleke, T.B., Si, B.C., 2006. Characterizing scale-dependent spatial relationships between soil properties using multifractal techniques. Geoderma 134 (3–4), 440–452.

stability of surface soil moisture in a semi-arid watershed. Remote Sens. Environ. 112 (2), 304–313. Dunne, T., Black, R.D., 1970. Partial area contributions to storm runoff in a small NewEngland watershed. Water Resour. Res. 6 (5), 1296 (&). Entin, J.K., Robock, A., Vinnikov, K.Y., Hollinger, S.E., Liu, S.X., Namkhai, A., 2000. Temporal and spatial scales of observed soil moisture variations in the extratropics. J. Geophys. Res.-Atmos. 105 (D9), 11865–11877. Evertsz, C.J.G., Mandelbrot, B.B., 1992. Self-similarity of harmonic measure on DLA. Phys. A Stat. Mech. Appl. 185 (1–4), 77–86. Famiglietti, J.S., Rudnicki, J.W., Rodell, M., 1998. Variability in surface moisture content along a hillslope transect: Rattlesnake Hill, Texas. J. Hydrol. 210 (1), 259–281. Feder, J., 1988. Fractals. Springer US, New York. Gomez-Plaza, A., Alvarez-Rogel, J., Albaladejo, J., Castillo, V.M., 2000. Spatial patterns and temporal stability of soil moisture across a range of scales in a semi-arid environment. Hydrol. Process. 14 (7), 1261–1277. Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. Oxford University Press. Goovaerts, P., 1998. Geostatistical tools for characterizing the spatial variability of microbiological and physico-chemical soil properties. Biol. Fertil. Soils 27 (4), 315–334. Grassberger, P., 1983. Generalized dimensions of strange attractors. Phys. Lett. A 97 (6), 227–230. Groenendyk, D.G., Ferré, T.P.A., Thorp, K.R., Rice, A.K., 2015. Hydrologic-process-based soil texture classifications for improved visualization of landscape function. PLoS One 10 (6), e0131299. Halsey, T.C., Jensen, M.H., Kadanoff, L.P., Procaccia, I., Shraiman, B.I., 1987. Fractal measures and their singularities: the characterization of strange sets. Nucl. Phys. B, Proc. Suppl. 2, 501–511. Hentschel, H.G.E., Procaccia, I., 1983. The infinite number of generalized dimensions of fractals and strange attractors. Physica D 8 (3), 435–444. Hu, Z.L., Islam, S., Cheng, Y.Z., 1997. Statistical characterization of remotely sensed soil moisture images. Remote Sens. Environ. 61 (2), 310–318. Huel, D., 2000. Managing Saskatchewan Wetlands: A Landowner's Guide. Saskatchewan Wetland Conservation Corporation, Regina, SK. Hupet, F., Vanclooster, M., 2002. Intraseasonal dynamics of soil moisture variability within a small agricultural maize cropped field. J. Hydrol. 261 (1–4), 86–101. Ji, W.J., Lin, M., Biswas, A., Si, B.C., Chau, H.W., Cresswell, H.P., 2016. Fractal behavior of soil water storage at multiple depths. Nonlinear Process. Geophys. 23 (4), 269–284. Journel, A.G., Huijbregts, C.J., 2003. Mining Geostatistics. Blackburn Press. Kachanoski, R.G., Dejong, E., 1988. Scale dependence and the temporal persistence of spatial patterns of soil-water storage. Water Resour. Res. 24 (1), 85–91. Kim, G., Barros, A.P., 2002. Downscaling of remotely sensed soil moisture with a modified fractal interpolation method using contraction mapping and ancillary data. Remote Sens. Environ. 83 (3), 400–413. Koster, R.D., Dirmeyer, P.A., Guo, Z.C., Bonan, G., Chan, E., Cox, P., Gordon, C.T., Kanae, S., Kowalczyk, E., Lawrence, D., Liu, P., Lu, C.H., Malyshev, S., McAvaney, B., Mitchell, K., Mocko, D., Oki, T., Oleson, K., Pitman, A., Sud, Y.C., Taylor, C.M., Verseghy, D., Vasic, R., Xue, Y.K., Yamada, T., Team, G., 2004. Regions of strong coupling between soil moisture and precipitation. Science 305 (5687), 1138–1140. Kravchenko, A.N., Bullock, D.G., Boast, C.W., 2000. Joint multifractal analysis of crop yield and terrain slope. Agron. J. 92 (6), 1279–1290. Li, Y., Li, M., Horton, R., 2011. Single and joint multifractal analysis of soil particle size distributions. Pedosphere 21 (1), 75–83. Mandelbrot, B.B., 1982. The Fractal Geometry of Nature. W.H. Freeman and Company, San Francisco.

23