Published in Soil Sci. Soc. Am. J. 68:577-587 (2004).
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
DIVISION S-6SOIL & WATER MANAGEMENT & CONSERVATION
Scale-Dependent Relationship between Wheat Yield and Topographic Indices
A Wavelet Approach
Bing Cheng Si* and
Richard E. Farrell
Dep. of Soil Science, Univ. of Saskatchewan, Saskatoon, Canada, S7N 5A8
* Corresponding author (Bing.Si{at}usask.ca).
 |
ABSTRACT
|
|---|
Topography can have a significant influence on crop yield, thus a better understanding of the effects of topographical parameters on crop yield is importantespecially for site-specific soil management. The objective of this study was to determine whether topographical indices developed for hydrological studies could be used as indicators of crop yield using a wavelet approach. The effects of soil curvature, upslope length, and a wetness index on wheat grain yields in a hummocky terrain were investigated along a transect on a clay loam Black soil (Blaine Lake Association) in Saskatchewan, Canada. A wavelet approach was used to elucidate the processes underlying the relationships between crop yield and topographical parameters. Wheat grain yields had significant correlations with upslope length (R2 = 0.60) and the wetness index (R2 = 0.46), whereas soil surface curvature explained only 15% of the variations in grain yield. Wavelet analyses revealed that significant variations occurred at scales <160 m for wheat grain yield, 280 m for upslope lengths, and 110 m for the wetness index. The cross wavelet analysis indicated that significant covariance existed at scales <180 m between wheat grain yield and upslope lengths and at scales <140 m between wheat yield and the wetness index, at a 99% confidence level.
 |
INTRODUCTION
|
|---|
THE SPATIAL VARIABILITY associated with crop yield is an integrated reflection of spatial variability in soil properties and factors affecting crop growth (i.e., light, temperature, and others). In rolling landscapes, topography is the most important factor controlling soil water redistribution, organic matter, nutrients, soil textural composition, and other soil properties that affect plant growth within a field (Zebarth and de Jong, 1989; Pennock et al., 1994; Moulin et al., 1994). It also affects variations in soil and canopy temperature and humidity. Hence, topographic parameters can be regarded as composite parameters that reflect the combined influence of various yield-affecting factors and their interactions. Consequently, an understanding of how soil variability affects the spatial variability of crop yields is a necessary prerequisite to developing proper soil management and precision farming techniques.
A valuable advantage of topographical data is that they are easy to obtain and relatively time-invariant compared with the measurement of more dynamic soil properties. Past studies have revealed a clear association between crop yield and topographical parameters (e.g., slope percentage, soil surface curvature, and elevation). For example, in years when precipitation is low to average, the influence of soil surface curvature on crop yield is reflected in higher yields in concave areas of the landscape (i.e., depressions) (Sinai et al., 1981; Timlin et al., 1998; Pennock et al., 2001; Kravchenko and Bullock, 2000). Conversely, water logging in the depressions may result in decreased crop yields during wet years (Timlin et al., 1998; Kravchenko et al., 2000). As well, elevation, slope, and aspect have been shown to exert a significant influence on wheat yield (Yang et al., 1998) and soybean quality (Kravchenko and Bullock, 2002). Nevertheless, associations between topography and crop yield are often weak in different fields. For example, Yang et al. (1998) found that selected topographic parameters could explain only 13 to 35% of the yield variability in wheat fields. Likewise, Miller et al. (1988) found no correlation between slope percentage and wheat yield.
In semi-arid or arid regions, soil water is the major limiting factor for crop production and processes that control the soil water distribution also control crop production. Water accumulation and runoff processes are largely determined by landscape configuration (convex vs. concave). Sinai et al. (1981) showed that in arid regions where there was significant unsaturated lateral soil water movement, soil water content was highly correlated with soil surface curvature. Zebarth and de Jong (1989) found that runoff during major precipitation events as well as the redistribution of snow and water from snowmelt resulted in soil water contents being greater in depressions than in knolls. That said, however, the spatial and temporal variability associated with soil water give rise to a level of complexity that prevents the estimation of soil water content based on a single factor, such as slope. Beven and Kirkby (1979) proposed a physically based composite wetness index to describe the topographic controls on soil water content by assuming that subsurface groundwater flow dominates soil water distribution and that any given point in the catchment is connected to its upslope area. This index has been used successfully to predict soil water deficit and surface saturation zones (Grayson and Blöschl, 2000; Bardossy and Lehmann, 1998). Gomez-Plaza et al. (2001) modified the procedure by including slope aspect to reflect another important factor in semiarid areasnamely, insolation. However, the revision of the Beven and Kirkby wetness index proposed by Gomez-Plaza et al. (2001) is not physically based and will not be used in this study. Another parameter is the contributing area or upslope length. Because snow redistribution and snowmelt runoff is the major processes for water redistribution in the rolling landscapes in semi-arid zone, upslope length may be a suitable wetness index in semi-arid zone. To date, however, there have been no reports of attempts to use the wetness indices such as upslope length, wetness index of Beven and Kirkby (1979)to predict the spatial variability associated with crop yields.
Topographical parameters, crop yield, and the relationships that exist between these variables exhibit both localized (Stevenson et al., 2001) and scale-dependent features (Miller et al., 1988). To characterize these features, researchers have employed an array of methods including: geostatisics (Miller et al., 1988), spectral analysis (Kachanoski et al., 1985), state-space methods (Wendroth et al., 1992), and multifractal approaches (Kravchenko et al., 2000). Whereas these methods allow complex relationships to be analyzed in either the spatial or the spectral domain, they cannot accommodate both domains simultaneously. Thus, because it is the totality of scale-dependent (spectral domain) and localized (spatial domain) features that determine a spatial series, new tools are needed to reveal both scale and localized information on yield and topographical parameters. Wavelet analysis represents one such tool that has been applied to problems in geophysics (Kumar and Foufoula-Goergiou, 1997), hydrology (Labat et al., 2000), and soil science (Lark and Webster, 1999, 2001; Si, 2003). Unlike Fourier analysis, which yields an average amplitude and phase for each harmonic in a data set, the wavelet transform produces an instantaneous estimate or local value for the amplitude and phase of each harmonic. This allows detailed study of nonstationary spatial or time-dependent signal characteristics. To the best of our knowledge, applications of the wavelet approach to the analysis of crop yield data have not been reported. The objective of this study was to examine the relationship between crop yield and selected topographic parametersincluding upslope length and the wetness index of Beven and Kirkby (1979)using wavelet analysis. We selected upslope length as a wetness index because upslope length is conceptually appealing for semi-arid zone where snowmelt runoff and snow redistribution are the water redistribution processes. We selected the wetness index of Beven and Kirkby (1979) because it is physically based and well accepted.
 |
THEORY
|
|---|
Wavelets
Information about the scales of the processes and spatial variability is important for understanding these processes in an agricultural field. A simple plot of a spatial data series shows all the variation in the data, that is, has the best spatial resolution, but shows nothing about the interrelationships between different length scales (Lindsay et al., 1996). Fourier transforms and the associated spectral analysis separate spatial variances at different scales, but in so doing, lose spatial information (spatial variability). For a spatial series collected in an agricultural field, it may have a trend (increasing clay content from one end of the field to another), or abrupt changes (i.e., hot spots in N2O emission). These trends and abrupt changes may indicate important soil, topographic, and ecological processes that affect crop production. However, traditional approaches are designed for stationary spatial series and new method such as the wavelet analysis may fill the gap.
Comprehensive reviews on the application of wavelet analysis can be found in Farge (1992) and Kumar and Foufoula-Georgiou (1997). A comparison between Fourier analysis and wavelet analysis is given in Kumar and Foufoula-Georgiou (1997) and Si (2003). In this paper, we only present basics regarding to wavelet analysis and help the readers familiarize with the interpretation of wavelet analysis. Wavelet is a small wave. A small wave grows and decays essentially in a limited distance. Sine function and cosine function are big waves because they keep on oscillating up and down on a plot of sin(x) or cos(x) versus x from
and
. To examine abrupt changes in space and frequency, the wavelet function must have a sufficiently rapid decrease around the origin of space, and an equivalently rapid decrease of Fourier transform of
(x) around the origin of the frequencies (localized in both space and frequency). In addition, this function must have a zero mean (Farge, 1992; Kumar and Foufoula-Georgiou, 1997). Mathematically, two basic properties (usually referred to as admissibility conditions) can be expressed as:
 | [1] |
and the square of
(x) integrates to unity:
 | [2] |
The second condition guarantees
(x) not to be zero everywhere and decrease to zero rapidly and first condition makes
(x) have equal negative and positive parts, such that the two parts cancel each other. There are numerous functions satisfying these conditions. Examples of such functions include the Harr, Mexican hat, and Morley wavelets (Fig. 1)
. In the following, we use Harr wavelet to illustrate the essence of continuous wavelet transform.

View larger version (30K):
[in this window]
[in a new window]
|
Fig. 1. The Harr wavelet and the Mexican hat wavelets for different scales and translation. (a) The Harr wavelet for scale = 1 (solid line) and 0.5 (dashed line); (b) the Mexican hat wavelets for scale = 1 (solid line) and 0.5 (dashed line); (c) the Harr wavelets for location = 0 (solid line) and 2 (dashed line); and (d) the Mexican hat wavelets for location = 0 (solid line) and 4 (dashed line).
|
|
Continuous Wavelet Transform
Assume y(x) is a spatial series where y is the measurement and x is the spatial location. The integral wavelet transform is defined by
 | [3] |
where
 | [4] |
and the asterisk corresponds to the complex conjugate. The function
(x), which can be a real or complex function, is called a wavelet. The
parameter can be interpreted as a dilation (
> 1) or contraction (
< 1) factor of the wavelet function
(x), corresponding to different scales of observation. The parameter,
, can be interpreted as a temporal or spatial translation or shift of the function
(x). If
= 0, we are looking at features at start of the data series. If
= 5 m, we are looking at features 5 m away from the start of the data series (the size of the feature we are looking at, depends on the scale parameter
). Thus
allows the signal, y(x), to be studied locally around the location
.
The Harr wavelet can be expressed as:
 | [5] |
It is not difficult to show that the Harr wavelet satisfies the conditions stated above. Harr wavelet is essentially composed of two Box car probability distribution functions; one is positive and one is negative. Therefore, for
= 1 and
= 0, we can separate the Harr continuous wavelet transform into two parts: one part being the integral of the product between x = 0.5 to x = 0, and another parts being the integral of the product between x = 0 to x = 0.5. Because integral of the product of probability with a function gives the average of the function, the first part is exactly the negative average of y(x) between 0.5 and 0 and second part the average between x = 0 to x = 0.5. Beyond the interval between x = 0.5 and x = 0.5, the Harr wavelet is zero, thus, integral of the product of y(x) and
(x) in this interval is also zero. Therefore, the Harr wavelet transform of y(x) is the average of y(x) at interval between x = 0.5 and x = 0, minus the average of y(x) at interval between x = 0 and x = 0.5. In another word, the Harr wavelet transforms extract information about how much difference there is between the two unit scale average of y(x) bordering on location x (for this example x = 0). Similarly, for the Mexican hat wavelets, the wavelet coefficient, W(
,
), at a specific
and
, is the sum (integral) of the product of y(x) and
(x), and can be considered as the difference of a weighted average of y(x) by
(x) between B and C, minus the weighted average of y(x) by
(x) between A and B, and minus the weighted average of y(x) by
(x) between C and D (Fig. 1b, dashed line). Beyond AD,
(x) = 0, hence, the product is also zero. Therefore, the Mexican hat wavelet coefficient, W(
,
), is proportional to the difference between a weighted average on unit scale and an average of two weighted average surrounding it. For a large
, the intervals AB, BC, and CD are larger and thus the averages are taken over a larger scale; however, a smaller weight [contracted
(x) value] is assigned to y(x) in calculation of the averages. Apparently, for spatially varying data series, different scale factor
, will results in different wavelet coefficients (or the differences in averages) for the same shift factor
. Obtaining a large collection of wavelet coefficients as a function of scales, allows us to examine the scale-dependence of the wavelet coefficients. Because
(x) has nonzero values only between A and D, the wavelet coefficient W(
,
) is not influenced by y(x) beyond AD, thus is localized between AD. By varying
, we shift the center of the wavelet from one location to another and obtain a wavelet coefficient associated with this specific location and scale. So, for a specified scale
we have a large number of wavelet coefficients corresponding to various translation parameters
. By varying
and
, we get a collection of wavelet coefficients as a function of scale
and location
. Therefore, the wavelet domain of a one-dimensional function (a spatial series along a transect) is two dimensional: one dimension corresponding to scale and the other to space (translation).
The Harr wavelet transform removes any linear trend in the spatial series automatically, because it takes the difference between consecutive segments in the data series (like a first order differencing) and the difference will not change with locations for linear trends. The Mexican hat wavelet transform, however, removes polynomial trends in the data series (the Mexican hat wavelet is a mini-polynomial). For the same reason, the Harr wavelet is very sensitive to abrupt increases or decreases and gives highest wavelet coefficients at sides of a peak, while the Mexican hat is sensitive to a peak in the data series, and gives highest wavelet coefficient at the peak location.
Different wavelet functions have different frequency and spatial resolutions and should be chosen for different applications. One common problem associated with wavelet analysis is the edge effect because the two ends in the data series like abrupt changes, and exert influences on points away from the ends for wavelet transform, particularly for large scale
. We refer the regions in which edge effect become important as cone of influence. The edge effects contaminate information in wavelet coefficients as a function of scales and locations, not only at the end, but also at locations away from the end, effectively reducing the length of the data series. The edge effects corroborate the problems of insufficient data length for application of wavelet analysis in agricultural and soil sciences. It is important to use a wavelet that has least edge effects. Of the often-used wavelet functions, the Mexican hat has less oscillations and is narrower than most of other wavelets. As a result, it has a much smaller cone of influence and is thus less affected by edge effects. In addition, Mexican hat is real-valued and capture both the positive and negative oscillations of the spatial series as separated peaks in wavelet power. The Morlet wavelet (Farge, 1992; Si, 2003) is complex and contains more oscillations than the Mexican hat, and hence the wavelet power combines both positive and negative peaks into a single broad peak. The separation of positive and negative peaks is important for analyzing the covariance of two variables, particularly at small scales. The Harr wavelet is very simple, but is not continuous. It is mainly used for education purpose. Therefore, we chose the Mexican hat wavelet, which can be expressed as:
 | [6] |
Where m = 2 and
() is the Gamma function. Since our measurements are usually taken at discrete points, continuous wavelet transform cannot be implemented directly. However, the continuous wavelet transform (Eq. [1]) can be implemented through discrete fast Fourier transform for improved efficiency. To adapt to the discrete form, we make the following changes of variables in Eq. [1]:
 | [7] |
where k and n are the location indices, s is the scale factor and
is the sampling intervals. To satisfy Eq. [2], the change of variables leads to:
 | [8] |
The wavelet function is expressed as a function of n, not x. Taking the Fourier transform of Eq. [3] and utilizing convolution theorem, we have,
 | [9] |
where a variable with a hat represents the Fourier transform of the variable.
The inverse Fourier transform of Eq. [9] gives
 | [10] |
where s is the wavelet scale, n is time index,
 | [11] |
and 
is the Fourier transform of the wavelet function. For the Mexican hat, 
is (Torrence and Compo, 1998):
 | [12] |
where
() is the Gamma function.
Wavelet Power Spectrum
The Fourier power spectrum measures the contribution of variance at a specific scale to the total variance and can be constructed from the sum of squares of Fourier coefficients at a certain frequency. Likewise, we can derive the wavelet power spectrum from the local wavelet coefficients. Because the wavelet function
(x) can be complex, the wavelet coefficient Wn(s) is also complex. The coefficient can then be separated into the real part and imaginary part, or amplitude, |Wn(s)|, and phase, tan1{
[Wn(s)]/
[Wn(s)]}. Consequently, one can also define the wavelet power spectrum as |Wn(s)|2. As indicated by Torrence and Compo (1998), for real-valued wavelet functions such as the Mexican hat, the imaginary part is zero and the phase is undefined. The global wavelet spectrum is the average of local wavelet spectra over all locations and is given as:
 | [13] |
where N is the number of data points in the spatial series. Percival (1995) showed that the global wavelet spectrum provides an unbiased and consistent estimation of the true power spectrum of a spatial series. Thus, the global wavelet spectrum provides an unbiased estimation of the contribution of variance at a specific scale to the total variance in a spatial series and different wavelet functions result in the same global wavelet spectrum.
The significance test for the global wavelet spectrum can be performed against a Gaussian white or red noise. A white noise is that the spatial series are independent identically distributed normal random variables, with zero mean and equal variance. The Fourier spectra of white noise do not change with scale, indicating that all scales (or frequencies) are equally present. In engineering application, a white noise is usually treated as the background and the null hypothesis for statistical test. A red noise is the data series that are dependent of each other in a short sampling interval, with a univariate lag-1 autoregressive process (Torrence and Compo, 1998). Soil properties are much similar at two nearby locations, because there are always some mixings of mineral particles at two nearby locations (due to sedimentation, biological activities, etc.), resulting in smaller variability for small sampling intervals. For large sampling intervals, those mixings are less likely, resulting in large spatial variability. Therefore, soil properties measured in the field have an increasing variance with increases in scales, resembling a red noise. This red noise-like behavior of a soil property does not depend on what field one is working on. Unless there are other processes superimposed on this red-noise background, the measured soil properties should behave like a red noise. If the measured soil properties deviate from the red-noise background, there must be other processes altered the red-noise process. Therefore, we can treat the red noise-like behavior as the background or the null hypothesis for statistical test. Because most geophysical series and agricultural field and data series can be modeled as autoregressive processes (Wendroth et al., 1992; Kachanoski and de Jong, 1988), we choose the red-noise background, which has been recommended by Muller and MacDonald (2000) and Torrence and Compo (1998) for geophysical applications. The discrete Fourier power spectrum of a red noise with unit variance is (Torrence and Compo, 1998; Shumway and Stoffer, 2000):
 | [14] |
Where N is the number of locations, k = 0, 1,..., N/2 is the frequency index, and r is the first order autocorrelation coefficient (Shumway and Stoffer, 2000). Clearly, the power spectrum of a red noise decreases with increase in frequency or increases with increase in scale. Assuming a red noise background spectrum, the distribution for the global wavelet power spectrum for the Mexican hat wavelet (Torrence and Compo, 1998), is best described by Chi square distribution
2a,d
k/d with a degree of freedom d:
 | [15] |
where
is the significance level. Clearly, the degree of freedom d is a function of wavelet scale, because local wavelet coefficients are not independent of each other. As a matter of fact, the larger the scale s, the stronger the correlation among wavelet coefficients and the smaller the degree of freedom d.
For significance test of global wavelet power spectra of a variable, the lag-1 autocorrelation coefficient for the variable should be calculated. The power spectrum of the red noise was then calculated using Eq. [14] and corresponding power spectrum
2a,d
k/d at a significance level of
is calculated with d given by Eq. [15]. If the actual global wavelet spectrum of a variable at a specific frequency or scale is larger than
2a,d
k/d, the difference between the global wavelet spectra of a variable and a red noise at that scale is statistically significant.
Cross Wavelet Spectrum
Given the two spatial series, y and z, with wavelet transforms W2n
and WYn
, the cross wavelet spectrum WYzn
can be defined as
 | [16] |
The cross wavelet spectrum is generally complex, hence we can define the cross-wavelet power as |WYzn
|. The confidence levels for the cross-wavelet power can be derived from the square root of the product of two chi-square distributions (Jenkins and Watts, 1968) and are given by Torrence and Compo (1998). The global cross wavelet spectrum is the average of local cross wavelet spectra over all locations and is given as:
 | [17] |
The global cross-wavelet power provides a measure of contribution of covariance at a specific scale to the total covariance.
Yzn
of a red noise has the following probability distribution with the degree of freedom provided by Eq. [15]
 | [18] |
where K0() is the modified Bessel function of order zero. For a confidence level of 99% (p = 0.99), Zd(p) can be calculated from Eq. [18]. If the calculated global cross wavelet spectrum >Zd(p)
/d, then
Yzn
of two variables is significantly different from that of red noises at a confidence level of p, where
Xk and
Yk are the power spectra of a red noise with lag-1 autocorrelation coefficient equal to that of variable X and Y at frequency index k, respectively.
Calculation of Curvature and Wetness Indices
Three topographic parameters were used to define the wetness index. The curvature was calculated from the elevations using the following approximations (Sinai et al., 1981; Pennock et al., 1994; Shary et al., 2002)
 | [19] |
The derivatives in Eq. [19] were approximated as
 | [20] |
where Z is elevation, i represents the indices for the x coordinates (along the transect) and j represents the indices for the y coordinates (perpendicular to the transect).
The wetness index of Beven and Kirkby (1979) was calculated as
 | [21] |
where
is the contributing area per unit contour length and tanß is the local terrain slope of the landscape elements. Hereafter, we refer to
as the upslope length, which is calculated as the distance from a given point in the landscape to the highest elevation point along the local slope.
 |
MATERIALS AND METHODS
|
|---|
Site Description
The study was performed in a semiarid area at Alvena, SK, Canada (45 km north of Saskatoon). The long-term (1971 to 2001) annual precipitation averages 350 mm, mostly falling in the spring (snow) and summer (rainfall). The average annual temperature is 2.2°C and potential evapotranspiration reaches 624 mm yr1; thus, the water deficit can be as high as 274 mm yr1. For 2001, the total precipitation was 159 mm, only 45% of the long-term average.
The landscape is classified as hummocky and the dominant soil type is an Aridic Ustoll. Five parallel transects (612 m in length) were established across the landscape with a lateral distance of 10-m between transects. A laser theodolite was then used to measure elevations at 6-m intervals along each transect (103 locations per transect). In all, a total of 515 points of elevation were measured in an area of 3.06 ha (i.e., 612 by 50 m). These elevation measurements allow for precise calculation of slope percentage, curvature, and upslope length at any point along the transect. The difference between minimum and maximum elevation along the transect was about 7 m; the maximum terrain slope was approximately 7°, or 12%.
Spring wheat was seeded on 26 Apr. 2001. All necessary management procedures were implemented as needed during the growing season. Fertilizer was applied at a uniform rate of 50 kg N ha1. Growing season precipitation obtained from a nearby weather station was 85 mm, which is about half the long-term average of 168 mm. On 30 Aug. 2001, a 1.0 by 1.0 m area near each measurement point along the center transect was selected and hand harvested for determination of total aboveground biomass and grain yield. The yield was converted to kilograms per hectare (kg ha1).
Statistics of measured grain yield, curvature, upslope length, and the wetness index of the 103 sites along the center transect are summarized in Table 1. Undisturbed soil samples (0 to 5 cm) were collected using a ring. Two sets of subsamples were used for determination of soil texture by the hydrometer method (Gee and Bauder, 1979) and organic C content using method of Wang and Anderson (1998). Clay content ranged from 20 to 41%; soil texture varied from loam to clay loam.
View this table:
[in this window]
[in a new window]
|
Table 1. Statistics of topographical parameters and grain yield obtained at 103 locations along the central sampling transect.
|
|
The continuous wavelet analysis was used to examine the spatial variability in wetness index, upslope length, and grain yield. To facilitate comparison between local wavelet spectra, crop yield, upslope length, and the wetness index were standardized by subtracting their mean from measurements and then dividing the difference by their standard deviation. Wavelet transform was implemented with the fast Fourier transform. Since Fourier transform operates on infinite data series and we have data series with limited data length, Fourier transform automatically connect the beginning of the data series with its end (wrap-around). Effect of wrap-around embedded in Fourier transform can artificially alter wavelet coefficients at the end of the data series. At the same time, fast Fourier transform require the length of the data series to be a power of two. There were 103 measurements for wheat grain yield, wetness index, and upslope length along the transect. Therefore, we padded 154 zeros to the end of wheat yield, the wetness index, and upslope length, to eliminate the wrap-around effect and to make the length of data series to be 256 (=28). The covariance relationships between grain yield and wetness index and between grain yield and upslope length also were analyzed using the cross-wavelet spectrum. All analyses were performed using a program written in Mathcad 2000 (Mathsoft Inc., Cambridge, MA) and utilizing the built-in fast Fourier transform.
 |
RESULTS AND DISCUSSION
|
|---|
Figure 2
shows the measured elevation and grain yield, along with the calculated slope curvature, upslope length, and wetness index along the transect. There are two large depressions (centered at 90 and 430 m) and a few small depressions centered at 200, 300, and 570 m along the transect. The measured grain yield, upslope length, and the wetness index exhibit a similar trend: large values in the depressions and small values on the knolls. In addition, measured grain yield and the upslope length increase sharply near the two large depressions at 90 and 430 m. The wetness index also increases near these depressions, though the increase is not as prominent as that of grain yield and upslope length. Clearly, spatial variations in the wheat yield and upslope length are nonstationary, exhibiting localized features at the two large depressions.
There was a strong correlation between grain yield and upslope length along the transect (R2 = 0.60***, Table 1). This is not surprising, given that water is the limiting factor for crop production in semiarid regions. For spring wheat, the number of heads is determined in the first month of wheat growth. Hence, water stress in spring can significantly reduce the number of heads per unit area, resulting in a loss of the wheat yield. Long-term average May precipitation was 3.3 cm, which is much less than the long-term monthly average potential evaporation (10.3 cm) in May. May precipitation in 2001 is 2.1 cm. Therefore, soil water storage before seeding becomes the most important water supply to spring wheat in May in the prairies, Canada. Soil water storage at spring is largely determined by the amount of snow accumulated over the winter and the snowmelt water accumulated during the snowmelt period, which has been documented by Zebarth and de Jong (1989), Kachanoski and de Jong (1988), and Pennock et al. (2001). In general, spring soil water storage at a given point in the landscape is closely related to the area contributing snow and water to that point. Thus, in dry to normal years, the larger the upslope length is, the greater the soil water storage and grain yield are. Upslope length also was correlated to the organic C content of the surface layer (R2 = 0.30***). Although soil organic C undoubtedly contributed to the high correlation coefficient obtained between upslope length and grain yield, the variation in grain yield explained by organic C content alone was much smaller than that explained by upslope length.
There was only a mild correlation between grain yield and wetness index (R2 = 0.46***), indicating that the wetness index explained less of the total variation in grain yield than did upslope length. This reflects the fact that the wetness index is equal to the ratio of the upslope length to the local slope at a given point in the landscape (Eq. [21]) and, as such, reflects the steepness of the slope at that point. This weighting is appropriate in humid regions where the main mechanism for water accumulation in low-lying areas is steady state, subsurface groundwater flow and runoff caused by water in excess of saturation (Beven, 1997). In semi-arid regions, however, snow redistribution and snowmelt runoff are the main mechanisms for water accumulation in low-lying areas (Zebarth and de Jong, 1989; Hayashi et al., 2003) and water accumulation at a point in the landscape is less sensitive to steepness of the slope than the upslope length. Therefore, in semi-arid regions, upslope length is a better indicator of soil water storage than either slope percentage or the wetness index.
There was a weak correlation between soil surface curvature and grain yield (R2 = 0.15***). Although this correlation was considerably lower than that reported by Sinai et al. (1981), it was consistent with the correlation coefficients reported by others (Timlin et al., 1998). The low coefficient of determination between soil surface curvature and grain yield presumably reflects the fact that soil surface curvature is independent of the magnitude of the upslope length and is essentially a measure of whether water flow at a point is convergent or divergent.
Scales of variation in grain yield, upslope length, and the wetness index, were analyzed using the wavelet approach. As suggested by Si (2003), the initial step involved exploratory analysis of the local wavelet spectrum to identify whether any patterns existed in the data and to determine whether these patterns were repeated across the transect (a global event) or were restricted to only one, or a few, localized regions across the transect (localized events). The periodicity (or frequency) of the repeated pattern is referred to as the scale of variation. If the pattern is global in scale, the global wavelet power spectrum is then analyzed to determine the significance of the variation.
Visual inspection of the local wavelet spectrum for grain yield (Fig. 3a)
revealed three scales of variation across the transect: 0 to 60 m (majority of region with low variances, or white color), 60 to 180 m (region interlaced with low and moderately high variances), and 180 to 300 m (300 m being the maximum scale resolvable by wavelet analysis for this study). For the 0- to 60-m scale, the variances were relatively uniform and small, indicating small scale-variance in the wheat yield was small. For the 60- to 180-m scale, regions of high variance were centered along the transect at 30, 90, 160, 200, 260, 300, 380, 430, and 520 m. These regions correspond to the locations of knolls and depressions across the transect and indicate the presence of a global feature. For the scales above 180-m, there were only four plumes of high variance centered at 90, 240, 430, and 570 m, which correspond to the locations of the largest depressions and knolls (Fig. 2a). To examine the distribution of variances in grain yield as a function of scale and perform statistical test of the global features for significance, the global wavelet power spectrum (Fig. 3b) was obtained by integrating the local wavelet spectra across the transect. The contribution of variance at a scale to the total variance increased with increase in scale. However, the curve at scales >160 m flattened off, suggesting that these scales have approximately equal contribution to the total variance. Spatial variations at scales <160 m were significantly different (at a confidence level of 99%) from that of a red noise spectrum (Fig. 3b). That is, variability at a spatial scale <160 m was not random. Conversely, variations at scales between 160 to 180 m were not significantly different (at a confidence level of 99%) from that of the red noise spectrum (i.e., were random). This indicates that the spatial variability associated with wheat grain yields was controlled by knolls and depressions with a scale <160 m. Variability above 180 m is localized features and global wavelet power spectrum cannot be used to test for statistical significance, because these localized features were not cycled events (Si, 2003).

View larger version (58K):
[in this window]
[in a new window]
|
Fig. 3. Local (a) wavelet spectrum and (b) global wavelet power spectrum of grain yield. Gray scale is expressed as the natural logarithm of the local wavelet spectrum. The dashed line in (b) is the power spectra of a red noise at a confidence level of 99%.
|
|
The local wavelet spectrum for upslope length (Fig. 4a)
revealed patterns similar to those observed for grain yield (Fig. 3a). That is, the local wavelet spectrum of upslope length also exhibited three scales of variations at 0 to 60, 60 to 180, and >180 m. For the 60- to 180-m scale, regions of high variance were centered at 20, 80, 150, 210, 250, 300, 360, 430, 540, and 600 m across the transect. Again, these locations correspond to the knolls and depressions occurring along the transect. At scales >180 m, the local wavelet power spectrum for upslope length also was very similar to that of grain yield, which is clearly a consequence of localized depressions and knolls. Since significance test can only be performed for global features in global wavelet power spectrum, significance test of these localized features (with scale >180 m for this case) was not made. At scales >180 m, the global wavelet power spectrum indicated that spatial variability associated with the upslope length was significantly different (at a confidence level of 99%) from that of the red-noise spectrum (Fig. 4b).

View larger version (54K):
[in this window]
[in a new window]
|
Fig. 4. Local (a) wavelet spectrum and (b) global wavelet power spectrum of upslope length. Gray scale is expressed as the natural logarithm of the local wavelet spectrum. The dashed line in (b) is the power spectra of a red noise at a confidence level of 99%.
|
|
The wavelet spectrum for the wetness index (Fig. 5a)
also exhibited features similar to those for grain yield (Fig. 3a) and upslope length (Fig. 4a). The scales of variations for the wetness index, however, occurred at 0 to 30, 30 to 180, and 180 to 300 m. Moreover, regions of high variance in the local wavelet spectrum of the wetness index encompassed a much smaller area than that of either grain yield or upslope length. The global wavelet power spectrum (Fig. 5b) revealed that spatial variability associated with the wetness index followed a trend similar to that observed for grain yield (Fig. 3b) and upslope length (Fig. 4b), but was significantly different from that of the red noise spectrum for a scale <110 m, at a confidence level of 99%), indicating that there was a distinct spatial pattern at these scales to the wetness index. At scales between 110 and 180 m, the variances were not significantly different (at a confidence level of 99%) from that of the red-noise spectrum. Therefore, wetness index of Beven and Kirkby (1979) smeared the features in upslope length between scales 110 and 180 m.

View larger version (60K):
[in this window]
[in a new window]
|
Fig. 5. Local (a) wavelet spectrum and (b) global wavelet power spectrum of the wetness index of Beven and Kirkby (1979). Gray scale is expressed as the natural logarithm of the local wavelet spectrum. The dashed line in (b) is the power spectra of a red noise at a confidence level of 99%.
|
|
Wavelet analysis also was employed to determine where a strong covariance existed between grain yield and upslope length or the wetness index. Figure 6a
shows the cross wavelet spectrum for grain yield and upslope length. At scales from 60 to180 m, strong covariance between grain yield and upslope length was manifested in the local cross wavelet spectrum centered at 30, 90, 150, 210, 260, 300, 380, 430, and 510 m along the transect. At scales >180 m, strong covariances between upslope length and wheat grain yield in the local cross wavelet power spectrum also were observed to correspond to the locations of the largest knolls and depressions (i.e., at distances of 100, 230, 430, and 540 m). To examine the distribution of covariance between grain yield and upslope length as a function of scale and perform statistical test of the global features for significance, the global cross wavelet power spectrum (Fig. 6b)obtained by integrating the local cross wavelet spectra across the transectwas significantly different (at a confidence level of 99%) from that of the red-noise spectrum at scales <180 m. Again, feature with a scale >180 m is localized, thus, cannot be tested for statistical significance.

View larger version (56K):
[in this window]
[in a new window]
|
Fig. 6. Local (a) wavelet cross-spectrum and (b) global wavelet cross-spectrum of grain yield and upslope length. Gray scale is expressed as the natural logarithm of the local wavelet spectrum. The dashed line in (b) is the power spectra of a red noise at a confidence level of 99%.
|
|
Figure 7a
shows the local cross wavelet spectrum for grain yield and the wetness index. Whereas the local cross wavelet spectrum is similar to that shown in Fig. 6a for grain yield and upslope length, the spatial pattern is not nearly as distinct. Areas with strong covariance in Fig. 7a are much smaller than those in Fig. 6a, indicating weaker covariance between grain yield and the wetness index than between grain yield and upslope length. At scales <140 m, the global cross wavelet power spectrum (Fig. 7b) was significantly different (at a confidence level of 99%) from that of the red noises. At scales >140 m, the cross wavelet spectrum for grain yield and the wetness index was not significantly different (at a confidence level of 99%) from that of the red noises.

View larger version (59K):
[in this window]
[in a new window]
|
Fig. 7. Local (a) wavelet cross-spectrum and (b) global wavelet cross-spectrum of grain yield and the wetness index of Beven and Kirkby (1979). Gray scale is expressed as the natural logarithm of the local wavelet spectrum. The dashed line in (b) is the power spectra of a red noise at a confidence level of 99%.
|
|
Significant scales of spatial variation in wheat grain yield and upslope length ranged from 0 to 160 and 0 to 180 m, respectively. The significant scale of covariation between wheat yield and upslope length was at scales <180 m. Therefore, it was concluded that wheat yields reflected a majority of the significant scales of variation associated with upslope length. For the wetness index, the local wavelet spectrum showed most of the same features as the spectrum for upslope length, though the spatial pattern for the wetness index was less distinct. According to Eq. [21], the smaller the upslope length at a given location, the smaller the wetness index and, hence, the less water stored in the soil. Small-scale depressions usually have small to medium slope steepness [tan(ß)]. Thus, the wetness index is dominated by the upslope length. On the other hand, medium-scale (<180 m) depressions can have high slope steepness, and the slope steepness becomes important. As a result, the spatial pattern discernible at scales >140 m for upslope length (and grain yield) was distorted for the wetness index. This may explain the smaller coefficient of determination between wetness index of Beven and Kirkby (1979) and grain yield than that between the upslope length and grain yield.
An important conclusion that can be drawn from this study is that upslope length exerted a significant influence on grain yield at scales <180 m. For scale from 0 to 180 m, the contribution of variance to the total variance increases with increase in the scale, indicating that the greater the upslope length, the higher the yield (positively correlated). Features within the scales of 0 to 180 m, are the dominant landscape features and cycled along the transect (global features). Localized features such as large depressions and knolls (with a scale >180 m) in this study can substantially affect crop yield and can be revealed through the local wavelet power spectrum of upslope and grain yield and cross-wavelet power spectrum between upslope length and wheat grain yield. Combination of global features and localized features provides the complete picture on the scale-location information of hydrological and agronomical processes in the field. Differentiation between localized and global features has the potential to provide guidance in designing efficient field soil management practices (based on global features) and site-specific soil management (based on localized features). The implications of these results are that: (i) farming practices should also consider the effect of upslope length, not only the landscape locations (i.e., knolls, backslopes, footslopes, and depressions); and (ii) wavelet analyses are useful for revealing localized and as well as global landscape features that exert significant controls on crop yields.
 |
SUMMARY
|
|---|
Statistical and wavelet analyses were used to examine the relationships between crop yield, upslope length, and wetness index at a site in the semi-arid agricultural region of Saskatchewan, Canada. Whereas the correlations between crop yield and both upslope length and the wetness index were significant, more of the spatial variability associated with crop yield was explained by upslope length than by the wetness index. Wavelet analysis was conducted to determine the scales of the variability in crop yield, upslope length, and wetness index. The local wavelet spectrum for crop yield was quite similar to that of upslope lengthand more so than to that of the wetness index. The global wavelet power spectrum was significantly different from that of a red noise at scales <160 m for crop yield, and 0 to 180 m for upslope length. At the same time, the local wavelet spectrum of the wetness index was significantly different from that of a red noise only at scales <110 m. The cross wavelet spectrum (used to assess covariance) between crop yield and upslope length showed significant differences from that of a red noise at scales <180 m at a 99% confidence interval. The larger the upslope length, the higher the yield is. This may explain why the slope curvature did not exert a strong influence on crop yield for this study because curvature reflects the shape, not the size of a depression.
Whereas spectral analysis provides scale information, localized features are lost (Si, 2003). And, without identifying these localized features, it is difficult to fully interpret and understand the scale information (such as the global wavelet power spectrum at scales >180 m). In the present study, wavelet analysis was used to reveal both scale information and the effect of localized features (namely, the large depressions) on wheat yield. In conclusion, it is suggested that wavelet analysis is an appropriate method for studying relationships between topographic indices and crop yield.
 |
ACKNOWLEDGMENTS
|
|---|
The research is funded by the National Science and Engineering Research Council of Canada (NSERC). The senior author thanks Drs. de Jong and Pennock for useful comments. Technical help from Shawn Campbell is appreciated. We thank Bob Hilkewich family for allowing us unlimited access to their farm field. We are grateful for the constructive comments from two anonymous reviewers.
Received for publication January 19, 2003.
 |
REFERENCES
|
|---|
- Bardossy, A., and W. Lehmann. 1998. Spatial distribution of soil moisture in a small catchment. Part 1: Geostatistical analysis. J. Hydrol. (Amsterdam) 206:115.
- Beven, K.J. 1997. TOPMODEL: A critique. Hydrol. Processes 11:10691085.
- Beven, K.J., and M.J. Kirkby. 1979. A physically-based, variable contributed area model of basin hydrology. Hydrol. Sci. Bull. 24:4369.
- Farge, M. 1992. Wavelet transform and their applications to turbulence. Ann. Rev. Fluid Mech. 24:395457.[Web of Science]
- Gee, G.W., and J.W. Bauder. 1979. Particle size analysis by hydrometer: A simplified method for routine textural analysis and a sensitivity test of measurement parameters. Soil Sci. Soc. Am. J. 43:10041007.[Abstract/Free Full Text]
- Gomez-Plaza, A., M. Martinez-Mena, J. Albaladejo, and V.M. Castillo. 2001. Factors regulating spatial distribution of soil water content in small semiarid catchments. J. Hydrol. (Amsterdam) 253:211226.
- Grayson, R., and G. Blöschl. 2000. Spatial patterns in catchment hydrology: Observations and modeling. Cambridge University Press. Cambridge, UK.
- Hayashi, M., G. van der kamp, and R. Schmidt. 2003. Focused infiltration of snowmelt water in partially frozen soil under small depressions. J. Hydrol. 270:214229.
- Jenkins, G.M., and D.G. Watts. 1968. Spectral analysis and its applications. Holden-Day.
- Kachanoski, R.G., D.E. Rolston, and E. de Jong. 1985. Spatial and spectral relationships of soil properties and microtopography: I. Density and thickness of A horizon. Soil Sci. Soc. Am. J. 49:804812.[Abstract/Free Full Text]
- Kachanoski, R.G., and E. de Jong. 1988. Scale-dependence and temporal persistence of spatial patterns of soil water storage. Water Resour. Res. 24:8591.
- Kravchenko, A.N., and D.G. Bullock. 2000. Correlation of corn and soybean grain yield with topography and soil properties. Agron. J. 92:7583.[Abstract/Free Full Text]
- Kravchenko, A.N., D.G. Bullock, and C.W. Boast. 2000. Joint multifractal analysis of crop yield and terrain slope. Agron. J. 92:12791290.[Abstract/Free Full Text]
- Kravchenko, A.N., and D.G. Bullock. 2002. Spatial variability of soybean quality data as a function of field topography: I. Spatial data analysis. Crop Sci. 42:804815.[Abstract/Free Full Text]
- Kumar, P., and E. Foufoula-Georgiou. 1993. A multicomponent decomposition of spatial rainfall fields: 1. Segregation of large- and small-scale features using wavelet transforms. Water Resour. Res. 29:25152532.
- Kumar, P., and E. Foufoula-Georgiou. 1997. Wavelet analysis for geophysical applications. Rev. Geophys. 35:385412.
- Labat, D., R. Ababou, and A. Mangin. 2000. Rainfall-runoff relations for karstic springs. Part II: Continuous wavelet and discrete orthogonal multiresolution analyses. J. Hydrol. (Amsterdam) 238:149178.
- Lark, R.M., and R. Webster. 1999. Analysis and elucidation of soil variation using wavelets. Eur. J. Soil Sci. 50:185208.
- Lark, R.M., and R. Webster. 2001. Changes in variance and correlation of soil properties with scale and location: Analysis using an adapted maximal overlap discrete wavelet transform. Eur. J. Soil Sci. 52:547562.
- Lindsay, R.W., D.B. Perceval, and D.A. Rothrock. 1996. The discrete wavelet transform and the scale analysis of the surface properties of sea ice. IEEE Trans. Geosci. Remote Sens. 34:771787.
- Miller, M.P., M.J. Singer, and D.R. Nielsen. 1988. Spatial variability of wheat yield and soil properties on complex hills. Soil Sci. Soc. Am. J. 52:11331141.[Abstract/Free Full Text]
- Muller, R.A., and G.J. MacDonald. 2000. Ice ages and astronomical causes: Data, spectral analysis and mechanisms. Springer-Verlag. New York.
- Moulin, A.P., D.W. Anderson, and M. Mellinger. 1994. Spatial variability of wheat yield, soil properties and erosion in hummocky terrain. Can. J. Soil Sci. 74:219228.
- Pennock, D.J., D.W. Anderson, and E. de Jong. 1994. Landscape-scale changes in indicators of soil quality due to cultivation in Saskatchewan, Canada. Geoderma. 64:119.
- Pennock, D.J., F. Walley, M. Solohub, B. Si, and G. Hnatowich. 2001. Topographically controlled yield response of Canola to nitrogen fertilizer. Soil Sci. Soc. Am. J. 65:18381845.[Abstract/Free Full Text]
- Percival, D.P. 1995. On estimation of the wavelet variance. Biometrika 82:619623.[Abstract/Free Full Text]
- Shary, P.A., L.S. Sharaya, and A.V. Mitusov. 2002. Fundamental quantitative methods of land surface analysis. Geoderma 107:132.
- Shumway, R.H., and D.S. Stoffer. 2000. Time series analysis and its applications. Springer-Verlag. New York.
- Si, B.C. 2003. Scale and location dependent soil hydraulic properties in nonlevel landscapes. p. 163178. In Y. Pachepski et al. (ed.). Scaling in soil physics. CRC Press, Boca Raton, FL.
- Sinai, G., D. Zaslavsky, and P. Golany. 1981. The effect of soil surface curvature on moisture and yieldBeer Sheba observation. Soil Sci. 132:367375.
- Stevenson, F.C., J.D. Knight, O. Wendroth, C. van Kessel, and D.R. Nielsen. 2001. A comparison of two methods to predict the landscape-scale variation of crop yield. Soil Tillage Res. 58:163181.
- Timlin, D., Y. Pachepsky, V.A. Snyder, and R.B. Bryant. 1998. Spatial and temporal variability of corn yield on a hillslope. Soil Sci. Soc. Am. J. 62:764773.[Abstract/Free Full Text]
- Torrence, T., and G.P. Compo. 1998. A practical guide to wavelet analysis. Bull. Am. Meteorol. Soc. 79:6178.
- Wang, D., and D.W. Anderson. 1998. Direct measurement of organic carbon content in soils by Leco CR-12 carbon analyzer. Commun. Soil Sci. Plant Anal. 29:1521.
- Wendroth, O., A.M. Alomran, C. Kirda, K. Richardt, and D.R. Nielsen. 1992. State-space approach to spatial variability of crop yield. Soil Sci. Soc. Am. J. 56:801807.[Abstract/Free Full Text]
- Yang, C, C.L. Peterson, G.J. Shropshire, and T. Otawa. 1998. Spatial variability of field topography and wheat yield in the Palouse region of the Pacific Northwest. ASAE 41:1727.
- Zebarth, B.J., and E. de Jong. 1989. Water flow in a hummocky landscape in Central Saskatchewan, Canada, I. Distribution of water and soils. J. Hydrol. (Amsterdam) 107:309327.
This article has been cited by other articles:

|
 |

|
 |
 
Q. Shu, Z. Liu, and B. Si
Characterizing Scale- and Location-Dependent Correlation of Water Retention Parameters with Soil Physical Properties Using Wavelet Techniques
J. Environ. Qual.,
October 23, 2008;
37(6):
2284 - 2292.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
X. Huang, L. Wang, L. Yang, and A. N. Kravchenko
Management Effects on Relationships of Crop Yields with Topography Represented by Wetness Index and Precipitation
Agron. J.,
September 8, 2008;
100(5):
1463 - 1471.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
C. L. Williams, M. Liebman, J. W. Edwards, D. E. James, J. W. Singer, R. Arritt, and D. Herzmann
Patterns of Regional Yield Stability in Association with Regional Environmental Characteristics
Crop Sci.,
July 1, 2008;
48(4):
1545 - 1559.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
B. C. Si
Spatial Scaling Analyses of Soil Physical Properties: A Review of Spectral and Wavelet Methods
Vadose Zone J.,
May 27, 2008;
7(2):
547 - 562.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
K. Parasuraman, A. Elshorbagy, and B. C. Si
Estimating Saturated Hydraulic Conductivity In Spatially Variable Fields Using Neural Network Ensembles
Soil Sci. Soc. Am. J.,
September 20, 2006;
70(6):
1851 - 1859.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
T. T. Yates, B. C. Si, R. E. Farrell, and D. J. Pennock
Wavelet Spectra of Nitrous Oxide Emission from Hummocky Terrain during Spring Snowmelt
Soil Sci. Soc. Am. J.,
May 23, 2006;
70(4):
1110 - 1120.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
A. N. Kravchenko, G. P. Robertson, K. D. Thelen, and R. R. Harwood
Management, Topographical, and Weather Effects on Spatial Variability of Crop Grain Yields
Agron. J.,
March 1, 2005;
97(2):
514 - 523.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
T. B. Zeleke and B. C. Si
Scaling Properties of Topographic Indices and Crop Yield: Multifractal and Joint Multifractal Approaches
Agron. J.,
July 1, 2004;
96(4):
1082 - 1090.
[Abstract]
[Full Text]
[PDF]
|
 |
|