Scale Dependence and Time Stability of Nonstationary Soil Water Storage in a Hummocky Landscape Using Global Wavelet Coherency

1. Introduction

Soil water is one of the most important limiting factors for semi-arid agricultural production and a key element in environmental health. It controls a large number of surface and sub-surface hydrological processes that are critical in understanding a broad variety of natural processes (geomorphological, climatic, ecological) acting over a range of spatio-temporal scales (Entin et al. , 2000). Knowledge on the behavior of soil water storage and its spatio-temporal distribution provides essential information on various hydrologic, climatic, and general circulation models (Beven, 2001; Western et al. , 2002), weather prediction, evapotranspiration and runoff (Famigleitti and Wood, 1995), precipitation (Koster et al. , 2004) and atmospheric variability (Delworth and Manabe, 1993).

The distribution of soil water in the landscape is the response of a number of highly heterogeneous factors and processes acting in different intensities over a variety of scales (Goovaerts, 1998; Entin et al. , 2000). The heterogeneity in factors and processes make the spatial distribution of soil water highly heterogeneous in space and time and create a challenge in hydrology (Quinn, 2004). Therefore, a large number of samples are needed in order to characterize the field averaged soil water with certain level of precision. However, if a field or watershed is repeatedly surveyed for soil water, some sites or points are consistently wetter or consistently drier than the field average. Vachaud et al. (1985) were the first to examine the similarity of the spatial pattern in soil water storage over time and termed this phenomenon time stability. The time stability is defined as a time invariant association between spatial location and classical statistical measures of soil water, most often the mean (Grayson and Western, 1998). Authors used the Spearman’s rank correlation to explain the similarity in the overall spatial patterns between two measurement series and the cumulative probability function of relative mean differences to examine the rank similarity of the individual locations over time [Vachaud et al. , 1985]. Various authors have used this concept to examine the similarity between the spatial patterns of soil water storage over a range of investigated area, sampling scheme, sampling depth, investigation period and land use (Kachanoski and de Jong, 1988; Grayson and Western, 1998; Hupet and Vanclooster, 2002; Tallon and Si, 2004; Martínez-Fernández and Ceballos, 2005; Starks et al. , 2006; Cosh et al. , 2008; Hu et al. , 2010b). However, information on the similarity between the spatial patterns of soil water within a season (intra-season), between seasons (inter-season), or within a season of different years (inter-annual) is not very common (Biswas and Si, 2011a).

Kachanoski and de Jong (1988) used the spatial coherency analysis to identify the similarity of the scales of the spatial patterns of soil water distribution over time and named the phenomena temporal persistence. Their study indicated loss of time stability at the scale < 40 m during the recharge period, which was attributed to topography. The spatial coherency analysis is based on the spectral analysis (Jenkins and Watts, 1968; Kachanoski and de Jong, 1988), which approximates the spatial data series by a sum of sine and cosine functions. Each function has an amplitude and a frequency or period. While the squared amplitude represents the variance contribution, the frequency component can be used to represent the spatial scale of ongoing processes (Webster, 1977; Shumway and Stoffer, 2000; Brillinger, 2001). The spectral analysis or frequency domain analysis is based on the assumption of second order stationarity (i. e. the mean and the variance of the series are finite and constant). However, more often than not, the soil spatial variation is nonstationary. Nonstationarity in the spatial distribution of soil water storage was also mentioned by Kachanoski and de Jong (1988). Nonstationarity restricts direct application of spatial coherency analysis to examine the similarity in the spatial patterns of soil water storage at different scales or scale-specific time stability, which calls for a new method.

Wavelet analysis can deal with localized features and thus nonstationarity by partitioning the spatial variations into locations and frequencies (Lark and Webster, 1999; Grinsted et al. , 2004; Si and Farrell, 2004, Yates et al., 2006; Biswas et al. , 2008), therefore providing an opportunity to study the spatial variation in soil water storage at multiple scales. While, the global wavelet analysis can deal with the scale specific variations, the global wavelet coherency analysis elucidates the scale specific correlation between any two spatial series. Therefore, the global wavelet coherency can be used to examine the similarity in the spatial patterns of soil water storage measured at two different times at multiple scales and study the scale-specific time stability. The objective of this study was to examine the scales of time stability of nonstationary soil water storage at different seasons in a hummocky landscape using the global wavelet coherency.

Figure 1.

Geographic location of the study site at St Denis National Wildlife Area within Prairie Pothole Region of North America along with the 3-dimensional and cross sectional view of the transect and different landform elements along the transect. CX indicates convex, CV indicates concave, CW indicates cultivated wetlands and UW indicates uncultivated wetlands

2. Theory

Wavelet analysis (Mallat, 1999) is used to divide a spatial series into different frequency components and study each component using a fully scalable window or wavelet. It calculates localized variations by shifting the standard function (mother wavelet) along the spatial series. The detail theory of the wavelet analysis is available elsewhere (Farge, 1992; Kumar and Foufoula-Georgiou, 1993, 1997; Torrence and Compo, 1998) and is beyond the scope of this chapter. There are different types of wavelet transform including discrete wavelet transform (DWT), continuous wavelet transform (CWT), wavelet packet transform (WPT), maximal overlap discrete wavelet transform (MODWT). These are suite of tools and can be used for certain purposes with some advantages and disadvantages. In this study, we use the continuous wavelet transform (CWT), where the wavelet coefficients at consecutive scales and locations can carry common information and provide a redundant representation of the signals information content and thus the detailed scale information (Farge, 1992; Lau and Weng, 1995; Keitt and Fischer, 2006; Furon et al., 2008). The detailed theory of the CWT can be found in various text books including Mallat (1998) and Chui (1992). Briefly, the CWT for a spatial series (Y_i) of length N (i= 1, 2, …,N) with an equal sampling interval of δx, can be defined as the convolution of Y_i with the scaled (s) and translated (x) wavelet (Torrence and Compo, 1998). Wavelet coefficients,$$W_i^Y(s)$$ can be calculated as

$$W_i^Y(s)=\sqrt{\frac{\delta x}{s}}{\displaystyle \sum _{j=1}^N{Y}_j}\psi \left[\left(j-i\right)\frac{\delta x}{s}\right]$$E1

where ψ[ ] denotes wavelet function. Out of many wavelet functions, the Morlet wavelet was used in this study because of enhanced spatial and frequency resolution. Morlet wavelet can be represented as (Torrence and Compo, 1998)

$$\psi (\eta )={\pi }^{-1/4}{e}^{i\varpi \eta -0.5{\eta }^2}$$E2

where, i is the complex number and equal to $$\sqrt{-1}$$, ω is the dimensionless frequency and η is the dimensionless space. The imaginary part conserved in the wavelet transform with Morlet wavelet can be used to identify the dominant orientation of variations in a random field. The energy associated with a scale and location can be measured from the magnitude of the wavelet coefficient. The wavelet power spectrum can be defined as $${\left|W_i^Y(s)\right|}^2$$, which is the space-frequency-energy representation of a spatial series. The global wavelet spectrum is the average of local wavelet spectra over all locations and is given by,

$${\overline{W}}^Y\left(s\right)=\frac{1}{N}{\displaystyle \sum _{n=0}^{N-1}{\left|W_n^Y\left(s\right)\right|}^2}$$E3

Similarly, the global wavelet spectra of another spatial series Z will be

$${\overline{W}}^Z\left(s\right)=\frac{1}{N}{\displaystyle \sum _{n=0}^{N-1}{\left|W_n^Z\left(s\right)\right|}^2}$$E4

The cross wavelet spectra between two spatial series Y and Z can be calculated as

$$W_i^{YZ}=W_i^Y{W}_i^Z{}^{\ast }$$E5

where, $$W_i^Z{}^{\ast }$$is the complex conjugate of $$W_i^Z$$. The global cross wavelet spectra can be calculated as

$${\overline{W}}^{YZ}\left(s\right)=\frac{1}{N}{\displaystyle \sum _{n=0}^{N-1}\left|W_n^{YZ}\left(s\right)\right|}$$E6

While the global wavelet cross spectra are similar to the covariances in the spatial domain, the global wavelet coherency spectra are similar to the coefficients of determination in the spatial domain for two variables. The global wavelet coherency spectra can be calculated as

$$R^2\left(s\right)=\frac{{\left|S\left({\overline{W}}^{YZ}\left(s\right)\right)\right|}^2}{S\left({\overline{W}}^Y\left(s\right)\right)S\left({\overline{W}}^Z\left(s\right)\right)}$$E7

where, S is the smoothing operator, $$S\left({\overline{W}}^{YZ}\left(s\right)\right)$$ is the smoothed global cross wavelet spectra of spatial series Y and Z, $$S\left({\overline{W}}^Y\left(s\right)\right)$$ and $$S\left({\overline{W}}^Z\left(s\right)\right)$$ are the smoothed global wavelet spectra of the spatial series Y and Z, respectively. In calculating wavelet coherency, it is necessary to smooth global cross wavelet spectra beforehand; otherwise, it will always be equal to 1 (Torrence and Compo, 1998; Maraun and Kurths, 2004). The coherency should be calculated on expected values. However, in most cases, there is only one realization of a spatial series, thus a coherency value has only one degree of freedom. By smoothing the coherency, one can overcome this problem and increase the degrees of freedom. In this study, we have used a boxcar window of size 5 (5 sample point average) to smooth the global wavelet and cross wavelet spectra.

Like the coefficient of determination, the global wavelet coherency spectra (R²) range from 0 to 1 and measure the correlation between two spatial series at each scale or within a particular frequency band. The closer the coherency values to one, the more similar the spatial patterns at a particular frequency or scale (= sampling interval / frequency).

The significance test for the wavelet coherency spectra can be carried out by calculating the confidence interval from an assumed theoretical distribution (Koopmans, 1974). However, the cutoff points for the test of hypothesis R² = 0 vs. R²> 0 can be conducted for s ≠ 0 from the F distribution (Koopmans, 1974):

$$R^2\left(s\right)=\frac{F_{2,N-2}\left(\alpha \right)}{n-1+F_{2,N-2}\left(\alpha \right)}$$E8

or

$$R^2\left(s\right)=1-{\left(1-\alpha \right)}^{\frac{1}{2m}}$$E9

where, α is the significance level and 2m + 1 is the width of the boxcar window. Therefore, m is the number of terms in each symmetrical half of the boxcar window. If the calculated coherency$${\widehat{R}}^2\left(s\right)$$ is greater than the theoretical value$$R^2\left(s\right)$$ at a particular scale (s), then the calculated coherency is significantly different from zero at the specified α. In this study, we have used m = 2, therefore the cutoff point at α = 0. 99 is 0. 684.

3. Materials and methods

A field experiment was conducted at St Denis National Wildlife Area, (52^°12′ N latitude, 106^°50′ W longitude), which is approximately 40 km east of Saskatoon, Saskatchewan, Canada (Fig. 1). A detailed information on the study site, soil water measurement and calibration of measurement instruments can be found in Biswas et al. (2012) and Biswas and Si (2011a, b). Briefly, the landscape of the study site is hummocky with a complex sequence of slopes (10 to 15%) extending from different sized rounded depression to complex knolls and knobs (Pennock, 2005) and is typical of the Prairie Pothole Region of North America (Fig. 1). The dominant soil type of the study site is Dark Brown Chernozem (Mollosiol in USDA soil taxonomy), which is developed from moderately fine to fine textured, calcareous, glacio-lacustrine deposits and modified glacial till (Saskatchewan Centre for Soil Research, 1989). The climate of the study area is semi-arid with the mean annual air temperature of 2^oC and the mean annual precipitation of 360 mm, of which 84 mm occurs as snow during winter (AES, 1997). The annual precipitation of the site during 2006, 2007, 2008, and 2009 were 489 mm, 366 mm, 331 mm, and 402 mm, respectively (Fig. 2). Year 2010 received 645 mm rainfall only during the spring and summer months (April to September), which is almost double the long-term average annual precipitation (Environment Canada, 2011).

Figure 2.

Monthly distribution of total precipitation in the year of 2006-09 along with the long term normal (90 year average)

Figure 3.

Site specific neutron probe calibration equation completed over three year time (2007-09). P indicates the ratio of the actual neutron count to the standard neutron count.

Soil water storage was measured along a transect of 576 m long with equally spaced 128 points (4. 5 m sampling interval). The transect was established over several knolls and seasonal depressions representing different landform cycles (Fig. 1). Topographic survey of the site was completed using Light Detection and Ranging (LiDAR) survey of the study area at 5 m ground resolution. Different landform elements were also identified as convex (CX), concave (CV), cultivated wetlands (CW) and uncultivated wetland (UW) (Fig. 1). The vegetation of the study site was mixed grass including Agropyronelongatum, Agropyronintermedium, Bromusbiebersteinii, Elymusdauricus, Festucarubra, Onobrychisviciifolia, Elymuscanadensis, Agropyrontrachycaulumand Medicago sativa, which was seeded in 2004 and allowed to grow every year. Surface 0-20 cm soil water was measured using time domain reflectometry (TDR) probe and a metallic cable tester (Model 1502B, Tektronix, Beaverton, OR, USA). A neutron probe (Model CPN 501 DR Depthprobe, CPN International Inc. , Martinez, CA, USA) was used to measure the soil water down to 140 cm at 20 cm vertical intervals. Soil cores at selected locations within 1 m around the neutron access tube were taken at different moisture conditions and the soil water content of each 10 cm interval were determined by gravimetric methods. The volumetric water content (gravimetric water content × bulk density) and the neutron counts were used to calibrate the neutron probe. The resulting calibration equation (θ_v = 0. 8523 P + 0. 0612 with n = 101 and r² = 0. 86, where P is the ratio of neutron count to standard neutron count) was used to convert the neutron probe count ratio to volumetric soil water at different depths and different locations (Fig. 3). Because neutron probe is prone to error for surface soil water measurements, the average soil water content at the surface 20-cm layer was measured using vertically installed time domain reflectometry probe and a metallic cable tester (Model 1502B, Tektronix, Beaverton, OR, USA). A standard calibration equation ($${\theta }_v=0.115\sqrt{k_a}-0.176$$, where $$k_a={\left[L_2/L\right]}^2$$ is the dielectric constant, L₂ is the distance between the arrival of signal reflected from the probe-to-soil interface and the signal reflected from the end of the probe curves (measured from waveform) and L is the length of the TDR probe) following Topp and Reynolds (1998) was used to derive the water content from the TDR recordings. Soil water content was measured 25 times at different seasons during a year over a five-year period (2007 – 2011). Based on the season, the measurements were divided into three groups: spring, summer, and fall. Though the analysis was completed for all the measurements in all the years, the space restriction in this chapter and the for demonstration purposes, only the result from 2008 and 2009 were used. The results from these years were very similar to the results from other years and can be generalized over the measurement period.

Wavelet analysis was completed using the MATLAB (The MathWorks Inc. ) code written by Torrence and Compo (1998) and is available online at http://paos. colorado. edu/research/wavelets/. The graphs were prepared in SigmaPlot (Systat Software Inc. ).

4. Results

High water storage was found at the locations of 100 to 140 m and 225 to 250 m from the origin of the transect, which were situated within depressions (Fig. 4). On contrary, the knolls stored less water. However, the difference of stored water in depressions and on knolls reduced from spring to fall within a year. For example, the range was 43. 42 cm for 20 April 2009 during spring and was 20. 81 cm for 27 October 2009 during fall season (Fig. 4). Apparently, the mean and variance of the first half of transect was quite different from that of the second half. It is evident that soil water along the transect is non-stationary.

Spearman’s rank correlation coefficients were used to examine the similarity of the overall spatial pattern. High rank correlation coefficients between any two-measurement series indicated time stability of overall spatial pattern of soil water storage. There was very strong intra-season time stability. For example, the rank correlation coefficient was 0. 98 between the measurement series on 2 May 2008 and 31 May 2008 (spring) and was 0. 99 between the measurement series on 23 August 2008 and 17 September 2008 (summer). Similarly, there was also strong inter-annual time stability. For example, the rank correlation coefficient was 0. 96 between the measurement series on 2 May 2008 (spring) and 7 May 2009 (spring) and was 0. 97 between the measurement series on 22 October 2008 (fall) and 27 October 2009 (fall). However, a relatively low rank correlation coefficient was observed between the measurement series from two different seasons. For example, the rank correlation coefficient was 0. 89 between the measurement series on 31 May 2008 (spring) and 23 August 2008 (summer), and 0. 85 between the measurement series on 31 May 2008 (spring) and 22 October 2008 (fall). However, the correlation coefficient was 0. 99 between the measurement series on 23 August 2008 (summer) and 22 October 2008 (fall) indicating strong similarity between summer and fall measurements. The correlation coefficients gradually decreased with the increase in time between measurements indicating the decrease in the degree of time stability over time. For example, the rank correlation coefficient was 0. 97 between the measurement series on 20 April 2009 and 7 May 2009, which gradually decreased to 0. 82 between the measurement series on 20 April 2009 and 27 October 2009.

Figure 4.

Spatial distribution of selected soil water storage series along the transect. The value in italics presents the average soil-water storage.

Figure 5.

Global wavelet coherency spectra of intra-season (spring, summer and fall of 2008) time stability. Values in the parentheses are the rank correlation coefficients.

Figure 6.

Global wavelet coherency spectra of inter-annual (spring, summer, and fall of 2008 and 2009) time stability. Values in the parentheses are the rank correlation coefficients.

Because of the nonstationarity of soil water along the transect, the scale specific similarity in the spatial pattern of soil water storage was examined using global wavelet coherency. There was strong intra-season time stability during summer or fall compared to spring (Fig. 5). Statistically significant strong coherency at all scales during summer and fall indicated that the spatial patterns present at different scales in summer were also present in fall. Similarly, large significant coherency between the measurement series of spring 2008 and spring 2009 or summer 2008 and summer 2009 indicated strong inter-annual time stability (Fig. 6). However, non-significant coherency at the scales < 20 m in the spring indicated the loss of intra-season time stability (Fig. 5). Similarly, there was loss of inter-annual time stability between springs of different years at the scales < 30 m (Fig. 6). However, the time stability was lost at scales < 65 m between the spring and summer and < 70 m between the spring and fall measurement series (inter-season) (Fig. 7). There was strong time stability between the summer and fall at all scales (Fig. 7). The minimum scale of statistically significant coherency was the lowest within a season and gradually increased with the increase in time between measurements (Fig. 8). For example, the minimum scale of significant coherency was 25 m between the measurement series on 20 April 2009 and 7 May 2009 and was 40 m between the measurement series on 20 April 2009 and 27 May 2009. The correlation between the measurement series on 20 April 2009 and 27 October 2009 was not significant at almost all scales except from 80 to 120 m (Fig. 8).

Figure 7.

Global wavelet coherency spectra of inter-season (spring-summer, spring-fall, and summer-fall of 2008) time stability. Values in the parentheses are the rank correlation coefficients.

5. Discussions

In our study area, depressions receive snowmelt runoff water from the surrounding uplands and store more water compared to knolls during the spring (Gray et al. , 1985; Winter and Rosenberry, 1995). In addition, uneven distribution of drifting snow in the landscape also contributes to the high water storage in depressions (Woo and Rowsell, 1993; Hayashi et al. , 1998; Lungal, 2009). Therefore, the alternate knolls and depressions along the transect created a spatial pattern in soil water storage inverse to the spatial pattern of elevation (Fig. 4). This spatial pattern with topography persisted through summer until fall (Fig. 4) as the depressions always stored more water than knolls. However, the variable demand of evapotranspiration reduced the difference in the maximum and minimum soil water storage (range) over time within a year. This phenomenon repeated every year.

Figure 8.

Global wavelet coherency spectra between the soil water measurement on 20 April 2009 and the measurements at different time within the year 2009. Values in the parentheses are the rank correlation coefficients.

At the absence of vegetation during spring, the soil water is lost mainly through surface evaporation, and to a lesser extent, the ground water interaction (Hayashi et al. , 1998; van der Kamp et al. , 2003). Evaporation may be higher in south-facing slope than in north-facing slope, but the difference in evaporation due to aspects may not be able to diminish the spatial patterns of soil water storage due to nonlocal controls (e. g. , macro-topography: knolls and depression) at large scales (Kachanoski and de Jong, 1988; Grayson et al. , 1997). However, the spatial patterns created from the micro-topography were not strong enough to dominate the differential evaporation created from the difference in local controls (e. g. , surface roughness, soil texture). Therefore, there was a loss of intra-season time stability at small scales during the spring. The loss of time stability at small scales was also observed between the spatial series from the spring of different years (Fig. 6). However, the influence of these small-scale processes was not strong enough to change the overall spatial pattern. Therefore, a strong intra-season and inter-annual time stability was observed in the overall spatial patterns (identified from Spearman’s rank correlation coefficients) irrespective of moisture conditions in different seasons (Martínez-Fernández and Ceballos, 2003).

However, with the establishment of vegetation, runoff from rainfall became the rare event and therefore the macro-topography was not able to reestablish the spatial patterns weakened by differential water uptake by vegetation (e. g. , more water uptake in depression than on knolls). Though macro-topography is still important, it is the interactions between vegetation and macro-topography that determined the spatial patterns of soil water. Because the interactions were relatively similar in later summer and early fall, the spatial patterns of soil water storage were very similar in summer and fall at all scales (Fig. 5). The intra-season and inter-annual time stability of the overall and the scale specific spatial patterns at all seasons contradicted the results of various authors who mentioned that the wet season showed strong time stability (Gómez-Plaza et al. , 2000; Qiu et al. , 2001; Hupet and Vanclooster, 2002), while others pointed out to the dry season (Robinson and Dean, 1993; Famiglietti et al. , 1998).

In addition, the processes controlling the spatial pattern of soil water storage may not change abruptly. For example, the vegetation growth is gradual and so is the evapotranspiration demand. In our study, large coherency between the measurement series on 20 April 2009 and 7 May 2009 indicated strong intra-season time stability. However, with the increase in time difference between measurements relatively weaker inter-season time stability was resulted (Fig. 8). This result contradicted the findings of Grayson et al. (1997), who indicated that the local controls (e. g. , vegetation, soil texture) are dominant in dry season (evapotranspiration > precipitation), while nonlocal controls (e. g. , topography) dominate redistribution of water during wet period or moisture surplus conditions (evapotranspiration < precipitation). However, the consistent coherency with reduced magnitude in our study indicated the change in the degree of the same control, but not a switch to different controls.

This study was different from Kachanoski and de Jong (1988), who used spatial coherency analysis to identify the scale specific similarity of the spatial patterns of soil water storage. The spatial coherency analysis assumed that the data series was stationary. However, authors identified the nonstationary nature of soil water storage (Kachanoski and de Jong, 1988) and divided the transect to create piecewise stationary series. Conversely, the global wavelet coherency analysis was able to deal with nonstationary soil water series. The wavelet analysis is well established to deal with nonstationarity. In addition, the conclusion of Kachanoski and de Jong (1988) was based on one-year measurement of soil water, which may not be universal as the precipitation variability over years may create a different experimental situation in different years. In this study, we have confirmed our conclusion based on five years of measurement of soil water storage.

Time stability is a result of multiple factors. Due to difference in the intensity of different factors, time stability of soil water and its scale dependence can be different. Instead of over-generalizing time stability, we classified time stability into intra-season, inter-season, and inter-annual time stability, because of the similar intra-season and inter-annual hydrological processes, but different inter-season hydrological processes. Therefore, the conclusion of this study may lead to improved prediction of soil water from reduced number of monitoring sites, thus allowing improved runoff and stream flow prediction in scarcely gauged basins.

6. Summary and conclusions

The similarity in the overall spatial pattern of soil water storage was first examined by Vachaud et al. (1985) and termed as the time stability of the spatial pattern. Kachanoski and de Jong (1988) extended the concept of time stability to the scale dependence of time stability using spatial coherency analysis. However, the stationarity assumption of the spatial coherency analysis restricts the use of this method for nonstationary spatial series. We have used global wavelet coherency analysis to examine the scale dependence of intra-season, inter-season and inter-annual time stability of nonstationary soil water spatial patterns.

There was strong intra-season time stability of the overall and scale specific spatial pattern. The time stability was lost at the scales < 20 m within the spring and < 30 m between the spring measurements from different years. However, strong time stability was present at all scales during the summer and fall, when the high evapotranspiration demand created similar spatial patterns. Similar processes in the summer and fall resulted strong inter-season time stability. However, not so similar processes in spring created weaker inter-season time stability between the spring and summer or the spring and fall. There was loss of time stability at the scales < 65 m and < 70 m between the spring and summer and the spring and fall, respectively. However, the change in the scales of time stability was not abrupt; rather it gradually decreased with the increase of time difference between measurements. The change in the similarity of the spatial patterns of soil water storage over time at different scales is an indicative of the change in the hydrological processes operating at those scales. Therefore, the analysis outcome can be used to identify the change in the sampling domain as controlled by the hydrological processes operating at different scales delivering the maximum information with minimum sampling effort.

Acknowledgments

The project was funded by a CSIRO Post Doctoral Fellowship and the University of Saskatchewan. Authors appreciate help from Amanda, Danny, Jason, Khizir and Henry in field data collection.