Dual scale trend analysis for evaluating climatic and anthropogenic effects on the vegetated land surface in Russia and Kazakhstan

We present a dual scale trend analysis for characterizing and comparing two contrasting areas of change in Russia and Kazakhstan that lie less than 800 km apart. We selected a global NASA MODIS (moderate resolution imaging spectroradiometer) product (MCD43C4 and MCD43A4) at a 0.05° (∼5.6 km) and 500 m spatial resolution and a 16-day temporal resolution from 2000 to 2008. We applied a refinement of the seasonal Kendall trend method to the normalized difference vegetation index (NDVI) image series at both scales. We only incorporated composites during the vegetative growing season which was delineated by start of season and end of season estimates based on analysis of normalized difference infrared index data. Trend patterns on two scales pointed to drought as the proximal cause of significant declines in NDVI in Kazakhstan. In contrast, the area of increasing NDVI trend in Russia was linked through the dual scale analysis with agricultural land cover change. The coarser scale analysis was relevant to atmospheric boundary layer processes, while the finer scale data revealed trends that were more relevant to human decision-making and regional economics.


Introduction
Widespread increases in plant growth across all northern latitudes were first reported in 1997 (Myneni et al 1997). Subsequent studies confirmed changes in spatial and temporal patterns of terrestrial vegetation, attributing these changes to warmer winters and springs (Nemani et al 2003). Direct human impacts on the land surface are especially evident in agricultural regions: 12% of the terrestrial surface is under active agricultural management (Ramankutty and Foley 1998, Cotton and Pielke 2007. Nevertheless, few global climate models incorporate the dynamic aspect of croplands and their sensitivities to anthropogenic influences (Osborne et al 2007). Global climate models instead typically represent agricultural regions by stable grassy vegetation that does not interact with the atmospheric boundary layer (Osborne et al 2007). Crops, however, display phenologies distinct from natural vegetation; the growing seasons within a region may be staggered in time; and crops may be irrigated (Vuichard et al 2008). In addition, crop establishment is generally fast and, in general, crops are rapidly removed at harvest (Bondeau et al 2007, Osborne et al 2007, Vuichard et al 2008. The areal extent and intensity of cultivation depend on local land use decisions that can be influenced strongly by distant institutions, policies, and markets (de Beurs and Henebry 2004a, Lerman et al 2004, Ioffe 2005, Lioubimtseva and Henebry 2009). The large amount of land surface used for crop production suggests that if agriculture and human geography are studied in isolation from climate change and variability, important feedbacks-both positive and negative-may be overlooked or underestimated (de Beurs and Henebry 2004a, Betts 2007, Osborne et al 2007. Thus, there is an urgent need to improve change attribution in studies of land surface dynamics. In particular, it is critical to consider direct human impacts on land surface dynamics originating from land use decisions. Earlier large scale vegetation studies were based on nine years of normalized difference vegetation index (NDVI) data derived from reflectance observations acquired by a series of advanced very high resolution radiometer (AVHRR) orbiting sensors (Myneni et al 1997). The NDVI exploits a spectral contrast between red and near infrared reflectance to indicate the presence of green vegetation (Tucker 1979). Changes in the vegetation type and/or density of the vegetated land surface can be expressed in terms of temporal trends in the NDVI retrieved from spaceborne sensors.
By the end of 2008, we acquired nine years of improved data (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) from NASA's moderate resolution imaging spectroradiometer (MODIS) sensor, which was first launched at the end of 1999 (on the Terra satellite) and again in 2002 in a different orbit (on the Aqua satellite). As a successor to the AVHRR, MODIS provides improved spatial and spectral resolution (Gallo et al 2005). A synoptic assessment of NDVI trends can reveal areas of change that merit closer attention because they indicate significant shifts in local water and carbon fluxes and in the surface energy balance (Kuemmerle et al 2008, Vuichard et al 2008, Henebry 2009).
Here we characterize and compare two contrasting areas of significant NDVI trends since 2000 that lie less than 800 km apart: positive trends in Russia versus negative trends in Kazakhstan.

Data processing
2.1. MODIS reflectance data: 0.05 • GCM grid and 10 • × 10 • tiles We selected a global NASA MODIS product (Terra + Aqua Nadir BRDF-Adjusted Reflectance data MCD43C4) at a 0.05 • (∼5.6 km) spatial resolution and a 16-day temporal resolution from the first available composite in 2000 through the last composite in 2008. We subset this global dataset to Eastern Eurasia, including Central Asia. In addition, we selected the same NASA MODIS product (Terra + Aqua Nadir BRDF-Adjusted Reflectance data; MCD43A4) at a 500 m spatial resolution and a 16-day temporal resolution from the first available composite in 2000 through the last composite in 2008. Higher resolution MODIS data (250-1000 m) are typically available as 10 • × 10 • lat/lon tiles with a sinusoidal projection. The globe exists of 36 tiles in horizontal direction and 18 tiles in vertical direction. These tiles are labeled in horizontal and vertical direction starting at tile h0v0. Here we have used twelve tiles: h20-23, v2-4. The total area included in the MODIS tiles incorporates 17 280 000 km 2 (or 69 120 000 pixels). The reflectance data in these products have been adjusted using models of bidirectional reflectance distribution functions to simulate reflectance from a nadir view (Lucht et al 2000). Both MODIS datasets are delivered as 8-day rolling composites, based on 16 days of data.
For each composite and tile we calculated two vegetation indices: the NDVI and the normalized difference infrared index or NDII. It has long been known that water strongly absorbs in the short-wave infrared (SWIR, MODIS band 7 = 2105-2155 nm) portion of the electromagnetic spectrum. SWIR reflectance is therefore sensitive to the amount of water in the vegetation (Gao 1996, Jang et al 2006. SWIR reflectance is generally low for high leaf water content, and increases with decreasing water content. The sensitivity of the SWIR channel to water has led to the development of a number of vegetation indices that are responsive to vegetation stress and vegetation removal based on SWIR and near infrared (NIR) reflectance (Jang et al 2006, Gu et al 2007. NDWI (Gao 1996) was developed from hyperspectral data as the difference between NIR reflectance and a SWIR band centered around 1240 nm. NDWI has been adapted for Landsat using the SWIR bands 5 (centered on 1650 nm rather than 1240 nm) or 7 at 2220 nm (e.g., Jang et al 2006, Gu et al 2007, Wang et al 2007. However, the index has been previously proposed as NDII using Landsat band 5 or the SPOT SWIR band (Hardisky et al 1983, Hunt andRock 1989) and as an index using Landsat band 7 (Hunt and Rock 1989). Higher NDII values are associated with higher water content in the leaves. The index is generally also high when there is snow on the land surface and drops during snow melt. In case of vegetation growth the index increases. We calculate the two indices as follows: NDVI = (NIR − Red)/(NIR + Red) and where NIR is MODIS band 2, Red is MODIS band 1, and SWIR is MODIS band 7. After calculating the two indices for each composite, we temporally resampled the 8-day product to 16 days by averaging consecutive 8-day composites. The final product consists of 12 MODIS tiles with 23 NDVI and NDII composites for each year (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) and 20 composites for 2000. A pixel time series was omitted from the analysis when either (1) it lacked more than 20% of the data or (2) it exhibited low NDVI seasonality, specifically, an average NDVI < 0.10 and a coefficient of variation of seasonal NDVI < 5%. These criteria filtered out deserts, inland water bodies, and cloudy and/or hazy areas; thus, excluding areas lacking either sufficient data or seasonality to produce reliable trend estimates. A total of 13.8% (table 1) of the land surface was excluded from analysis by this filtering.

MODIS land cover data
There are many different land cover datasets available. We selected the MODIS land cover data product (MOD12) from the Boston University group (http://duckwater.bu.edu/lc/ mod12q1.html) because it is created from the same satellite sensors as the NDVI dataset used for trend analysis. This product is based on MODIS data from 2001, is delivered globally at a 1 km spatial resolution, and includes several land cover classification schemes. We selected the IGBP global land cover classification scheme because it provides the largest number of classes. We subset the data for our study regions in Russia and Kazakhstan.

Landsat data
We selected global land survey (GLS) 2005 data in Russia and in Kazakhstan. These Landsat images correspond to some of the areas of significant positive and negative trends, respectively, in our analysis. The Russian study area is located at the most southern edge of Kirov oblast, going into Mari El oblast to the south. This region has been identified with depressed agriculture due either to suboptimal climatic conditions or to substantial drops in population (Ioffe et al 2006). Depressed agricultural areas were identified as regions with very low productivity both with respect to crops and livestock, such as low average grain yield, low number of cattle, and low average milk yield per cow (Ioffe et al 2004). The Kazakh study area is located at the northwestern border between Kazakhstan and Russia and consists mostly of grasslands with some areas of marginal agriculture.

NDII start of season (SOS) and end of season (EOS) determination
We have previously defined the phrase 'land surface phenology' (LSP) to refer to the spatio-temporal development of the vegetated land surface as revealed by satellite sensors (de Beurs and Henebry 2004a). Due to the spatial resolution of satellite sensors, LSP deals with mixtures of land covers and thus is distinct from the traditional notion of a species-centric phenology (Friedl et al 2006). LSP metrics are primarily based on image time series of vegetation indices (VI) from optical sensors. If a VI derived from satellite observations is to be used to monitor the duration of vegetation activity, it is desirable to compare the satellite retrieved phenological estimates with data observed at ground level. However, a principal disadvantage of phenological observation by satellite imagery is the complexity of validation of the data by ground observations that usually measure something quite different (Schwartz et al 2002, Fisher et al 2006, Liang and Schwartz 2009). As a result, it is often unclear what the LSP metrics actually track. For example, in high latitude biomes, the greatest temporal increase in the normalized difference vegetation index (NDVI), which some methods use to indicate 'start of the season' (SOS), is often due to snow melt (Reed et al 1994, Delbart et al 2005. The end of the season (EOS) metric that is based on a decline in the VI, on the other hand, can be confused with an extended period of cloudiness that yields lower NDVI, instead of actual senescence. Since the relationship between satellite measures of LSP and phenological events of particular species or lifeforms is ambiguous, a diversity of satellite measures and methods has arisen. The range of methods can be divided into four main categories: threshold, derivative, smoothing algorithms, and model fit.
In a recent paper we provided an intercomparison of a range of methods to determine SOS (White et al 2009). The midpoint NDVI method by White et al (1997) provided SOS estimates which were most consistent with available field observations in North America. However, this method might still be negatively affected by snow cover. Thus, we applied an approach based on the NDII index to determine the nominal start and end of the growing season (Delbart et al 2005(Delbart et al , 2006. Note that this method was not part of the intercomparison project, because it is based on data incorporating the SWIR band which was not available in the selected NOAA AVHRR dataset. Delbart et al (2005Delbart et al ( , 2006) used a related index, the normalized difference water index (NDWI), based on the near infrared (780-890 nm) and short-wave infrared (1580-1750 nm) bands of SPOT Vegetation. The underlying idea is that the index (NDII or NDWI) first decreases with snow melt, and then increases during initial vegetation growth. Thus, this approach is especially powerful to distinguish actual increases in vegetation from the loss of snow cover, which also appears as an increase in NDVI. Specifically, SOS is determined as the last date at which NDII is lower than the NDII minimum plus 20% of the total increase. EOS is determined as the date for which NDII has decreased by 20% of the fall amplitude. We note that this method has not been tested by Delbart et al (2006) in agricultural areas. However, we evaluated a range of threshold percentages and determined that 20% still provided the most consistent result for our study region (data not shown).

Trend analysis
Using simple linear regression to estimate a trend from a time series is a practice widespread in the remote sensing literature. However, there are four basic assumption that can affect the validity of trends summarized in a linear regression equation: (1) all ordinate values (i.e., mapped on the y-axis) should be mutually independent; the residuals should be (2) random with (3) zero mean; and (4) the variance of the residuals is independent of time. NDVI time series typically violate the first assumption (there is usually high positive autocorrelation between consecutive observations) and when multiple sensors are involved, the second and the fourth assumption are often violated as well.
We have previously discussed the seasonal Kendall (aka seasonal Mann-Kendall) trend test corrected for autocorrelation as an alternative to trend analysis by simple linear regression Henebry 2004b, 2005b).
The original Mann-Kendall (MK) trend test is nonparametric and is calculated by summing the number of times a particular observation has a higher value than any of the previous observations (Hirsch et al 1982). If the value of a particular composite is higher than a previous composite, one is added to the test statistic, if the values are equal, nothing is added and if the value is lower than a previous composite, one is subtracted from the test statistic. The seasonal Kendall (SK) trend test for image time series first calculates the MK statistic for each composite separately. The SK statistic for the complete time series consists of the sum of the MK statistics for all composites.
The test statistic is asymptotically normal with a zero mean. The variance is defined as the sum of variances for every 'season' (in our case, every composite) plus the sum of the covariances for every combination of 'seasons'. The White areas in the map have a retrieved length of season less than 50 days. Tan areas are masked out due to a lack of data or a lack of seasonality. (Bottom) Vegetation trends from 2000 to 2008 revealed by NASA MODIS sensors at 0.05 • spatial resolution. Areas outlined in orange and green indicate highly significant ( p 0.01) negative and positive trends, respectively. Areas in tan were excluded from analysis. Areas in shades of gray did not exhibit highly significant trends. Underlying grayscale image is the average NDVI over the study period, darker areas indicate lower average NDVI. The yellow boxes indicate the areas in Russia and Kazakhstan selected for finer scale analysis. autocorrelation correction is applied to the calculation of the covariance which is corrected for autocorrelation between consecutive seasons or composites (Hirsch and Slack 1984, Hess et al 2001, de Beurs and Henebry 2004b: where n is the number of years in the time series, r gh Spearman's rank correlation coefficient for composites g and h, and i and j the years. We calculate the SK trend test separately for every pixel for all composites during the growing season, where the growing season is determined as the average start of season and the average end of season (as determined by NDII) between 2001 and 2007. Thus, the number of composites incorporated in the trend analysis for each pixel is allowed to vary spatially, so northern locations generally incorporate fewer composites during the growing season than southern locations. (We did not incorporate the years 2000 and 2008 in the SOS/EOS analysis, because the first images in 2000 did not start in January and the year 2008 was still lacking a few winter images at the time of our analyses.) We applied the SK trend test to the temperature station data from March through September for 2000-2006 and to the gridded GPCC precipitation data to assess temperature and precipitation changes in the study areas over the past seven to eight years.

Results
We have divided the results according to the scale of the imagery. We discuss the results of the trend analyses first at a coarser spatial resolution (0.05 • ) and then at a grain size (500 m) more than 100× finer. To guard against spurious trends, we present only highly significant trends having a p-value of less than 0.01. Figure 2  (bottom) displays the highly significant trends ( p < 0.01) in the vegetated land surface since 2000. We found significant negative trends in 6.1% of the study region (2.1 × 10 6 km 2 ), with extensive change in Russia east of the Urals and across northern Kazakhstan (figure 2, table 1). Substantial percentages of negative trends can also be found in Estonia, Finland, Kyrgyzstan and Mongolia, although these countries only account for a small percentage of the study area land surface. Highly significant positive trends appeared in only 2.4% of the study area (0.82 × 10 6 km 2 ) with change areas in far northeastern Russia and China (table 1). Substantial percentages of positive trends can also be found in Armenia, Azerbaijan, Georgia and North Korea, although these countries only account for a small percentage of the study area land surface (table 1). Table 2 gives the most important land cover classes for each change category. For example, in Kazakhstan we found a positive trend in less than 1% of the country (7073 km 2 , table 1). These positive trends were predominantly located in the barren class (48.2% or 3409 km 2 , table 2) and in the open shrublands (18.1% or 1280 km 2 ). On the other hand we found a negative trend in 17.9% of the Kazakh land surface (487 649 km 2 , table 1). Of these negative trends 66.6% (or 324 774 km 2 ) were located in grasslands and 24.4% (or 118 986 km 2 ) were located in croplands. Thus, it appears that a small part of the barren and open shrubland areas in Kazakhstan saw an increase in NDVI, while large parts of the typically semi-arid grasslands and croplands saw significant declines in NDVI. A similar pattern can be found, for example, in Mongolia, where significant negative trends occur mainly in grasslands, while positive trends occur in barren classes. Table 2 shows that while certain countries reveal increases and decreases in the same classes (e.g., Azerbaijan has both positive and negative trends in croplands), other countries, such as Kazakhstan and Mongolia as discussed above, reveal very different classes in the positive and negative change categories.

Fine scale analysis
Attribution is the key challenge in any change analysis. Previous coarse resolution change maps have shown anomalies attributed to disturbance events (Potter et al 2003), temperature and precipitation changes, and correlations with large scale climate processes (Nemani et al 2003, Potter et al 2008. This change map based on the past nine years presents a mixture of positive and negative trends resulting from both direct and indirect impacts of climate variability and change. Human land use decisions also drive many of the observed changes. While several areas present interesting change patterns, here we will focus only on the significant negative changes found in Kazakhstan and the much smaller area of significant positive changes found just north of Kazakhstan in the agricultural areas of Russia. While the fine scale analysis has been performed for a large study area (figure 3), we discuss in detail the observed changes that occur within two Landsat scenes. We selected these two regions because: (1) summer Landsat images were available for both regions around the same time (August 2006 and August 2007) and (2) there was prior evidence of change. Figure 3 shows the results of the trend analysis at 500 m resolution. The trend pattern between the 0.05 • and 500 m resolution is very consistent with widespread negative changes in the grassland and shrubland areas of Kazakhstan and smaller areas of positive trends in Russia. The tan areas in the drier deserts of Kazakhstan were masked due to a lack of seasonality in the data. Figure 4 zooms into the selected areas in Russia and Kazakhstan, respectively, to illustrate the retrieved growing season length at 500 m resolution. On average for the selected  study regions, the growing season length is only 10 days shorter in Russia than Kazakhstan. Thus, the difference in the number of 16-day composites incorporated for the two regions is less than one. As a result, we do not anticipate the trend results between the two regions to be significantly affected due to the difference in incorporated composites.
To understand the nature of the observed trends, we superimposed the 500 m trend results on both the 1 km MODIS  land cover map and a 30 m Landsat image. Figure 5 gives the results for the selected Russian study region. Both the land cover map and the Landsat image reveal that the observed NDVI increases are predominantly confined to the cropland parts of the region, while the adjacent forested areas do not exhibit change. Figure 6 shows the corresponding results for the Kazakh study area. Declines are more widespread and span across both grassland and cropland classes, but the declines appear to be focused in the southern (grassland) parts of the area.
To investigate whether the observed changes were predominantly driven by weather, we analyzed the precipitation and temperature station data within the two study areas. In Russia, we find no significant temperature trend between 2000 and 2008 and no significant precipitation trend (based on the GPCC data) between 2000 and 2007 (table 3). In Kazakhstan, we do not find a significant temperature trend between 2000 and 2006 in any of the three stations. All three Kazakh stations, however, reveal a negative precipitation trend, although the trend is only significant for one of the stations at the 0.10 level.

Kazakhstan
For the total study region, the area of significant NDVI declines is almost three times as large as the area with significant NDVI increases. The largest NDVI declines appear to be located in Central Asia (Kazakhstan, Kyrgyzstan, Tajikistan, Turkmenistan, and Uzbekistan) where 13.6% of the land surface experienced NDVI declines over the past few years and just about 0.5% experienced NDVI increases. The regions with negative trends are very widespread and span multiple land cover classes (figure 6). These changes appear to be related to droughty conditions, especially in 2004, that occurred in the spring wheat regions and arid grasslands of Kazakhstan (Lindeman 2005) and in other Central Asian countries as well. FAO production statistics for Kazakhstan (figure 7) support this conclusion. The decline in wheat area cultivated, between 1992 and 2000, illustrates the socio-economic impact of institutional changes following the collapse of the Soviet Union. In wheat areas the fallow period extended indefinitely giving rise to weeds and grasses that were not heavily grazed due to a significant concomitant decline in livestock (de Beurs and Henebry 2004a). The wheat decline in the period before 2000 was predominantly a socio-economic effect (de Beurs and Henebry 2004a). Figure 7 reveals an increase in wheat area in the period between 2000 and 2007, indicating agricultural recovery in Kazakhstan. However, the yield and production variability reveal the effects of weather and yield and production were down sharply in 2003, 2004, and 2005. Negative precipitation trends were also confirmed by our station (gauge) analysis and the GPCC data (table 3). Others also confirmed these trends by analyses of recent station and GPCC precipitation data for Kazakhstan (Akhmadiyeva and Groisman 2008).

Russia
Agricultural reform has been one of the most important change processes in European Russia that has been unfolding since the formal collapse of the Soviet Union at the end of 1991 (Ioffe 2005, Ioffe et al 2006. Widespread land abandonment is perhaps the most striking consequence of the reform, visible even in synoptic imagery (de Beurs and Henebry 2004a, 2005b, Kuemmerle et al 2008. Land abandonment in Russia is not occurring randomly or, in fact, unexpectedly; it is preceded by persistently low crop yields (Ioffe et al 2006). While Russia is often treated as a unified whole, especially with respect to agricultural reform, substantial regional diversity exists (Wegren et al 2004). Our trend analysis indicates NDVI increases especially in the croplands classes and not in the other land cover classes, suggesting an anthropogenic rather than a climatic cause. Nevertheless, the analysis presented here cannot determine unequivocally whether this increase is due to land abandonment or to an increase in agricultural production. However, the selected study area in Russia is located in a region characterized by a depressed agricultural sector and emigration from rural to urban environment (Ioffe et al 2004(Ioffe et al , 2006. Furthermore, we have shown previously in Kazakhstan that agricultural land abandonment can lead to increased NDVI due to a succession of weeds, grasses and forbs, and eventually shrubs and trees that tend to green-up before annual crops would be planted Henebry 2004a, 2005a).

Comparing trends detected at two scales
We have revealed significant positive and negative NDVI trends in Russia and Kazakhstan at two vastly different spatial resolutions (5.6 km versus 500 m). While the general pattern of trends is comparable between the two scales, it is important to note that the trends detected at the coarser resolution could result from aggregates of patches of significant trends detectable only at finer resolutions. The study area in Russia demonstrates this scaling phenomenon especially well: the finer grained trends are associated with agricultural land cover-even specific fields or clusters of fields-while this patchy distribution of change is smeared across a larger area or even vanishes at the coarser resolution due to land cover mixtures within pixels. It is important to note that the rescaling of trend results is not straightforward due to (1) high spatial heterogeneity at the finer scale, (2) the nonlinearity of vegetation indices, and (3) the thresholding effect of specific significance levels. Furthermore, the coarser scale of the MODIS CMG (0.05 • or 5.6 km) is relevant to atmospheric boundary layer processes and, indeed, is at the effective resolution limit of current regional climate and mesoscale meteorological models. Thus, the smearing of the trends does not necessarily constitute a loss of information because there is some de minimis scale of landscape patchiness below which the boundary layer fails to respond to changes. The finer scale analysis reveals trends that are more relevant to human decision-making and regional economics (Ioffe et al 2004(Ioffe et al , 2006. Thus, even though the changes at the finer scale may disappear at the courser scale resolution as a result of mixed pixels, it would not be right to argue that these finer scale changes are not relevant or significant. The dual scale of the trend analysis enables a partitioning of change attribution that would be very difficult at a single scale.

Conclusions
Previously we argued for the SK trend test as a good alternative for often applied ordinary least squares trend analysis within sensor periods (de Beurs and Henebry 2004b). In this paper we have refined the method by determining trends only for the vegetative growing season as delineated by SOS and EOS determination based on NDII. Trend analysis is, however, primarily an exploratory tool that can highlight areas of interest, namely, those exhibiting significant change. Trend attribution remains a critical but challenging exercise. Our dual scale trend analysis, however, provides complementary views of the land surface and enables a partitioning of proximal causes. Here we use multiple lines of evidence to infer drought as the proximal cause of significant NDVI declines in Kazakhstan. The area of NDVI increase in Russia was initially indicated by the coarse scale trend analysis, but the finer scale analysis reveals the link with agricultural land cover change. This specific pattern of land cover change appears to be more related to a direct anthropogenic source of change than a climatic one. We point out that the underlying cause of some anthropogenic change could be climate change, e.g., people could abandon land due to severe droughts, resulting in a mixture of climatic and anthropogenic effects. However, our current conclusion is buttressed by other sources that indicate this area as a problematic with respect to agriculture (Ioffe et al 2006). Nevertheless, additional socio-economic analysis as well as time series of satellite data at even finer spatial resolution (e.g., 30 m Landsat data) will be required to address more fully the observed change dynamics.