Characterizing Cropland Phenology in Major Grain Production Areas of Russia , Ukraine , and Kazakhstan by the Synergistic Use of Passive Microwave and Visible to Near Infrared Data

We demonstrate the synergistic use of surface air temperature retrieved from AMSR-E (Advanced Microwave Scanning Radiometer on Earth observing satellite) and two vegetation indices (VIs) from the shorter wavelengths of MODIS (MODerate resolution Imaging Spectroradiometer) to characterize cropland phenology in the major grain production areas of Northern Eurasia from 2003–2010. We selected 49 AMSR-E pixels across Ukraine, Russia, and Kazakhstan, based on MODIS land cover percentage data. AMSR-E air temperature growing degree-days (GDD) captures the weekly, monthly, and seasonal oscillations, and well correlated with station GDD. A convex quadratic (CxQ) model that linked thermal time measured as growing degree-days to accumulated growing degree-days (AGDD) was fitted to each pixel’s time series yielding high coefficients of determination (0.88 ≤ r2 ≤ 0.98). Deviations of observed GDD from the CxQ model predicted GDD by site corresponded to peak VI for negative residuals (period of higher latent heat flux) and low VI at beginning and end of growing season for positive residuals (periods of higher sensible heat flux). Modeled thermal time to peak, i.e., AGDD at peak GDD, showed a strong inverse linear trend with respect to latitude with r2 of 0.92 for Russia and Kazakhstan and 0.81 for Ukraine. MODIS VIs tracked similar seasonal responses in time and space and were highly correlated across the growing season with r2 > 0.95. Sites at lower latitude (≤49◦N) that grow winter and spring grains showed either a bimodal growing season or a shorter unimodal winter growing season with substantial inter-annual variability, whereas sites at higher latitude (≥56◦N) where spring grains are cultivated exhibited shorter, unimodal growing seasons. Sites between these extremes exhibited longer unimodal growing seasons. At some sites there were shifts between unimodal and bimodal patterns over the study period. Regional heat waves that devastated grain production in 2007 in Ukraine and in 2010 in Russia and Kazakhstan appear clearly anomalous. Microwave based surface air temperature data holds great promise to extend to parts of the planet where the land surface is frequently obscured by clouds, smoke, or aerosols, and where routine meteorological observations are sparse or absent.


Introduction
About 10% of the earth's terrestrial surface (1.5 billion ha) is covered by crops, both irrigated and rainfed [1].In the face of increasing concerns about the food-water-energy nexus, a better understanding of global crop dynamics is urgently needed.Although grain production in Russia (RU), Ukraine (UA), and Kazakhstan (KZ) declined significantly following the collapse of the Soviet Union, it has begun to recover [2], and wheat production was largely unaffected by Russia's annexation of Crimea due to favorable weather [3][4][5].If the recent increases in crop production continue, this region may become a key source to address food security crises in the coming decades, despite ongoing challenges arising from land use change, regional conflicts, and climatic variability, extremes, and change.
The theoretical basis of using remotely sensed information to monitor crop growth and development was based on United States Department of Agriculture (USDA) research in the 1960s through the 1970s that characterized the optical properties of crop type and crop fields [6].The USDA Foreign Agricultural Service now routinely uses remote sensing to monitor conditions in key crop production areas across the globe to generate market intelligence, particularly in regions where information is scarce or unreliable [7].The normalized difference vegetation index (NDVI; [8]) is the most commonly used and a well-documented vegetation index.Remotely sensed vegetation indices (VIs) such as the NDVI have been widely utilized for agricultural mapping and monitoring [9][10][11][12][13].The United States Agency for International Development (USAID) funded and United States Geological Survey (USGS) implemented Famine Early Warning System Network (FEWS NET) uses NDVI-based measures of cropland activity as part of its integrated early warning system for food security and drought monitoring [14].FAO's Global Information and Early Warning System (GIEWS) uses NDVI data to detect vegetation health as a proxy for crop production.Despite the great advantages, the NDVI has some limitations including reduced sensitivity in denser vegetation, reduced detection of green vegetation in sparsely vegetated areas due to soil background, and reliance on just two spectral bands, the red and near infrared (NIR) bands [15][16][17].The enhanced vegetation index (EVI; [18]) was designed to overcome this loss of sensitivity over dense vegetation, adjust for soil background effects, and reduce atmospheric effects in the red band by including information from the blue portion of the spectrum [18].These vegetation indices-NDVI and EVI-are complementary for global vegetation studies and together improve detection of changes in surface vegetation and extraction of canopy biophysical variables [18,19].
Land surface phenology (LSP) deals with the timing of vegetated land surface dynamics as observed by satellite remote sensors, particularly at spatial resolutions and extents relevant to meteorological processes in the atmospheric boundary layer [20,21].LSP plays an important role in monitoring cropland dynamics.Comparative studies of cropland LSPs enable the discrimination between crop types [22][23][24][25], provide observations which can be assimilated into process-based crop models increasing model prediction accuracy [26,27], enable impact assessment and early warning for drought [28], and can aid forecasting of crop yields [29,30].
Phenological studies integrate climate-biosphere relationships because the timing of vegetation lifecycle events is influenced by temperature and precipitation [31][32][33][34].The mid-latitudes are thermally limited with respect to net primary production; temperature is the primary constraint for crop production [30,35,36].Air temperature is a key forcing for crop growth and development [37][38][39][40].The concept of heat units or thermal time was introduced by Réaumur in 1730 [41].A common metric of thermal time is the growing degree-day (GDD), which weights the passage of calendar time by the temperature deemed useful for plant growth [20,38,42].Use of GDD and accumulated GDD (AGDD) can improve phenological analyses [41,43,44], and crop yield models [45][46][47].
In most countries, observations from the network of weather stations are sparse in space and incomplete in time.Remote sensing offers multiple approaches to estimate temperature and moisture in a spatially comprehensive way at the cost of spatial precision and temporal resolution.Land surface temperature (LST) can be retrieved from moderate resolution (≤1 km) spatial resolution thermal infrared (TIR) sensors, such as MODIS (MODerate-resolution Imaging Spectroradiometer).Although LST is often very different than surface air temperature measured nominally at 2 m, some studies argue that it may be more relevant for modeling plant growth [48][49][50].However, TIR sensors cannot see through clouds to retrieve LST [20,51,52].Passive microwave radiometers, which operate at much coarser spatial resolution (25 km) relative to TIR sensors, can penetrate most clouds to retrieve air temperature, although not over frozen surfaces [52][53][54][55].Although 25 km may be considered very coarse relative to many kinds of land surface features and phenomena, it is still a much finer spatial resolution for air temperature data than exists across most of the global land surface [56,57].
Here we explore whether the synergistic use of vegetation indices from visible and near infrared (VNIR) reflectance and surface air temperature data retrieved through microwave radiometry can improve the characterization of cropland phenology at the mid-latitudes.The higher spatial but lower temporal resolution VNIR VIs are sensitive to green-up and brown-down dynamics of the growing season, and the coarser spatial but finer temporal resolution surface air temperature data provide the biometeorological tempo to improve modeling the land surface dynamics.In Section 2, we describe the data and modeling methods used.The results of our analyses appear in Section 3, including modeling the land surface seasonality of thermal time, latitudinal patterns in thermal time, and the identification of unimodal and bimodal seasonality.Section 4 discusses the identification of the period of canopy evapotranspiration detected from the thermal time residuals, and the detection of anomalous patterns arising from the impact of heat waves on crops.Section 5 concludes with remarks about the significance of this study and directions for further research.

Remote Sensing Data
We used products from two sensors: MODIS for the VNIR and TIR data and AMSR-E (Advanced Microwave Scanning Radiometer on EOS) for the microwave data.AMSR-E, onboard the NASA-EOS Aqua satellite since May 2002, is a multi-frequency microwave radiometer (6.925, 10.65, 18.7, 23.8, 36.5, 89.0 GHz) that detects faint microwave emissions from the earth's surface and atmosphere.The Numerical Terradynamic Simulation Group at the University of Montana has produced various geophysical parameters from AMSR-E twice daily data records (daytime, ~13:30, and nighttime, ~01:30) overpasses, and 25 km spatial resolution from mid-June 2002 to October 2011 (antenna failure).These include surface air temperatures (ta; ~2 m height), fractional coverage of open water over land (fw), vegetation canopy transmittance (tc) at three microwave frequencies (6.925 GHz, 10.65 GHz and 18.7 GHz), surface soil moisture (mv; ≤2 cm soil depth), and integrated water vapor content (V) of the total column [58].Only surface air temperature retrievals were used in this study [57].The average surface air temperature calculated from the 01:30/13:30 overpass pair is quite likely to differ from average air temperature calculated from the daily extrema.We do not consider this potential bias to be a problem for our analysis because the accumulated average daily temperature serves to provide the tempo for growth rather than establish temperature threshold that would trigger phenophase transitions.In other words, the discrepancies introduced by the potential bias are negligible with respect to the modeling.
MODIS is aboard the Terra and Aqua satellites launched by NASA in 1999 and 2002, respectively.MODIS has 36 spectral bands from which different groups of data products are processed at different spatial and temporal resolutions tuned for specific applications.For this study, we used MODIS collection 5 level 3 data products in a 0.05 degree (~5.6 km) Climate Modeling Grid (CMG).These products included 8-day composites of land surface temperature (LST) from both satellites (Terra: MOD11C2; Aqua: MYD11C2) and 16-day composites of Nadir BRDF-Adjusted Reflectance (NBAR), which combines observations from Terra and Aqua (MCD43C4; [59]).NBAR provides improved retrievals of surface reflectance through consistent normalization of multiple views of the surface to a nadir view using bidirectional reflectance distributions functions (BRDF) to model surface anisotropies [60,61].
To identify our specific cropland study areas, we used two datasets: (1) the International Geosphere Biosphere Programme (IGBP) land cover scheme in the MODIS land cover product at a spatial resolution of 0.05 • (MCD12C1; [59]); and (2) the USDA Foreign Agriculture Service (FAS) crop layers by oblast (for Russia: [62], for Kazakhstan: [63], for Ukraine: [64]).Details on the selection process of the study sites are presented in Section 2.2.

In Situ Data
Meteorological station data from National Oceanic and Atmospheric Administration (NOAA) historical climate data [14] was used to compare with the satellite retrievals.Among the 49 sites on which the study is conducted, we selected six meteorological stations within 25 km distance from the geographic centers of these study sites.From these six stations, a considerable amount of data was missing except for the year 2003 (Table 1).Missing data pose a problem for calculating accumulated growing degree-days.While there are more meteorological stations and data available within these countries, we restricted the analysis to what is freely available through WMO data-sharing agreements.In this study, we first selected our study cropland sites using criteria described in Section 2.2, and only then searched for meteorological stations within a reasonable distance from the selected cropland study sites.

Study Region
Our study region spans the major grain producing areas of Ukraine, southern Russia, and northern Kazakhstan.According to MODIS Land Cover Type 1 Percentage Product for the IGBP scheme (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010), the dominant classes in the study region are Grassland (IGBP class 10 at 34%), Cropland (IGBP class 12 at 27%), Crop-Natural Vegetation Mosaic (IGBP class 14 at 14%), and Mixed Forest (IGBP class 5 at 11%).The overall land cover classification accuracy for this product is about 75%, but the range in class-specific accuracies is large [61].In order to select specific AMSR-E study pixels with a high and stable percentage of cropland, we conducted a temporal stability analysis for each of the three dominant land cover classes over the 8 years study period by first calculating the maximum, minimum, and mean land cover percentage over the study period and then displaying the maximum, mean, and range in the red, green, and blue color planes, respectively [50].Yellow (high maximum and high mean and low range) identifies land cover that is stable over the study period (Table 2, Figure 1).Expanses of yellow in Figure 1 indicate where Crop-Natural Vegetation Mosaic (Figure 1a), croplands (Figure 1b), and grasslands (Figure 1c) continue to be the dominant land cover during the study period.1) arises from the false color composite of red, green, and blue color planes that display, respectively, the maximum percentage of LC class, the average percentage of LC class, and the range of percentages of LC class over the study period.Source: [50].For legend, refer to Table 2.
We selected the study sites within the stable core areas of cropland, but we also used the USDA Crop Explorer maps to identify wheat producing areas.We identified 49 AMSR-E pixels (study sites) in the three countries: 14 in Ukraine, 24 in Russia, and 11 in Kazakhstan (Figure 2, Table 3).Land cover percentage for each AMSR-E pixel study site was determined by aggregating the MCD12C1 Land cover percentage within each AMSR-E pixel.The collection of sites spans nearly 12° of latitude and 56° of longitude.The most extreme latitudinal sites (Cherkessk and Kazan', both in Russia) have a maximum day length difference of two hours at the summer and winter solstices.The site-to-site multi-year average cropland cover ranges between 46% (Volgograd, Russia) to 100% (a number of sites mostly in Ukraine), while the mean cropland cover for all sites over the study period was 88%.While the sites are dominated by crops, there is some degree of fragmentation or heterogeneity within the 25 km pixels.However, the Soviet era heritage of very large field sizes persists across the region and minimizes that fragmentation.The rest of the land cover proportion is mainly crop-natural vegetation mosaic (CNVM), followed by grassland.Factors that could affect the temporal variation in cropland cover include class indeterminacy due to site adjacency and spectral mixing, drought and recovery, soil freeze-thaw state, disturbance, and land cover change.For legend, refer to Table 2.
We selected the study sites within the stable core areas of cropland, but we also used the USDA Crop Explorer maps to identify wheat producing areas.We identified 49 AMSR-E pixels (study sites) in the three countries: 14 in Ukraine, 24 in Russia, and 11 in Kazakhstan (Figure 2, Table 3).Land cover percentage for each AMSR-E pixel study site was determined by aggregating the MCD12C1 Land cover percentage within each AMSR-E pixel.The collection of sites spans nearly 12 • of latitude and 56 • of longitude.The most extreme latitudinal sites (Cherkessk and Kazan', both in Russia) have a maximum day length difference of two hours at the summer and winter solstices.The site-to-site multi-year average cropland cover ranges between 46% (Volgograd, Russia) to 100% (a number of sites mostly in Ukraine), while the mean cropland cover for all sites over the study period was 88%.While the sites are dominated by crops, there is some degree of fragmentation or heterogeneity within the 25 km pixels.However, the Soviet era heritage of very large field sizes persists across the region and minimizes that fragmentation.The rest of the land cover proportion is mainly crop-natural vegetation mosaic (CNVM), followed by grassland.Factors that could affect the temporal variation in cropland cover include class indeterminacy due to site adjacency and spectral mixing, drought and recovery, soil freeze-thaw state, disturbance, and land cover change.

Site
Sites adjacent to land cover transition zones experience higher temporal variation in cropland cover (sites closer to the unstable peripheral areas with the magenta color in Figure 2 (see also Table 2)).For example, Saratov 1 in Russia (site 21) had the largest grassland encroachment with 69% and 31% mean cropland and grassland cover, respectively, over the eight years and, thus, had the largest cropland cover range (28%).Omsk 1 in Russia (site 46), and Kostanay 2 in Kazakhstan (site 36) are other examples of this phenomenon.

Methods
We used the AMSR-E and MODIS data from all eight years (2003-2010) for all the pixels.To avoid the frozen season and to maintain a consistent analysis period across all pixels, we restricted our focus to data from DOY 89-305, which ranges 30 March to 1 November for non-leap years and 29 March to 31 October for leap years.For the AMSR-E data, 8-day moving average filter was applied to the daily data to minimize data gaps due to orbit and swath width.Growing degree-day (GDD) is the daily thermal-time increment above a certain threshold (base temperature) for plant growth [20,41,45], and accumulated growing degree-days (AGDD) are the simple summation of heat units throughout the annual observation period [21,65,66]; that is, "the passage of days is weighted by the quantity of growing degrees occurring that day, with zero (but not negative) degrees being a permissible weight" [67].McMaster [68] concluded that the base temperature of 0 °C is a robust and sufficient threshold for wheat phenology, which is the main crop in our study area.GDDs were calculated for both AMSR-E air temperature (ta) data and MODIS land surface temperature (LST) data with a base temperature of 273.15K (=0 °C) as follows: where tASC and tDSC are the air temperature retrievals at the ascending and descending passes.GDDs were accumulated from DOY 89-305 to generate each year's AGDD trajectory.The GDDs calculated from MODIS LST 8-day composites were multiplied by eight to rescale them into the proper range for eight daily values, rather than a single value occurring in eight days.The daily AMSR-E GDDs were summed over eight days to align temporally them with the particular MODIS compositing periods.GDDs from both sensors were compared with available nearby meteorological station (stations within 25 km distance from the geographic centers of the AMSR-E pixel study sites) daily air temperature GDD to check for consistency using linear regression model.The station data is processed using similar methods as the AMSR-E daily air temperature data.We then test whether the slopes and/or intercepts from the linear regression model fit were significantly different by sensor.For this test, we used analysis of covariance ([ANCOVA; 69])-in which sensors (AMSR-E and MODIS) were the independent categorical effect variables, the meteorological station GDD was the covariate, and satellite GDD was the response variable.In this analysis, p ≤ 0.05 was considered significant.
Accumulated GDD (or AGDD) was calculated from GDD for both sensors.Sites adjacent to land cover transition zones experience higher temporal variation in cropland cover (sites closer to the unstable peripheral areas with the magenta color in Figure 2 (see also Table 2)).For example, Saratov 1 in Russia (site 21) had the largest grassland encroachment with 69% and 31% mean cropland and grassland cover, respectively, over the eight years and, thus, had the largest cropland cover range (28%).Omsk 1 in Russia (site 46), and Kostanay 2 in Kazakhstan (site 36) are other examples of this phenomenon.

Methods
We used the AMSR-E and MODIS data from all eight years (2003-2010) for all the pixels.To avoid the frozen season and to maintain a consistent analysis period across all pixels, we restricted our focus to data from DOY 89-305, which ranges 30 March to 1 November for non-leap years and 29 March to 31 October for leap years.For the AMSR-E data, 8-day moving average filter was applied to the daily data to minimize data gaps due to orbit and swath width.Growing degree-day (GDD) is the daily thermal-time increment above a certain threshold (base temperature) for plant growth [20,41,45], and accumulated growing degree-days (AGDD) are the simple summation of heat units throughout the annual observation period [21,65,66]; that is, "the passage of days is weighted by the quantity of growing degrees occurring that day, with zero (but not negative) degrees being a permissible weight" [67].McMaster [68] concluded that the base temperature of 0 • C is a robust and sufficient threshold for wheat phenology, which is the main crop in our study area.GDDs were calculated for both AMSR-E air temperature (ta) data and MODIS land surface temperature (LST) data with a base temperature of 273.15K (=0 • C) as follows: where t ASC and t DSC are the air temperature retrievals at the ascending and descending passes.GDDs were accumulated from DOY 89-305 to generate each year's AGDD trajectory.The GDDs calculated from MODIS LST 8-day composites were multiplied by eight to rescale them into the proper range for eight daily values, rather than a single value occurring in eight days.The daily AMSR-E GDDs were summed over eight days to align temporally them with the particular MODIS compositing periods.GDDs from both sensors were compared with available nearby meteorological station (stations within 25 km distance from the geographic centers of the AMSR-E pixel study sites) daily air temperature GDD to check for consistency using linear regression model.The station data is processed using similar methods as the AMSR-E daily air temperature data.We then test whether the slopes and/or intercepts from the linear regression model fit were significantly different by sensor.For this test, we used analysis of covariance (ANCOVA; [69])-in which sensors (AMSR-E and MODIS) were the independent categorical effect variables, the meteorological station GDD was the covariate, and satellite GDD was the response variable.In this analysis, p ≤ 0.05 was considered significant.
Accumulated GDD (or AGDD) was calculated from GDD for both sensors.
where GDD t is daily temperature increment of growing degree-days at time t.During summer in the Northern hemisphere, the maximum and minimum air temperatures typically occur late afternoon and predawn, respectively, while data from AMSR-E and MODIS Aqua (N.B.: AMSR-E is onboard Aqua satellite) are recorded twice daily (ascending-daytime, ~13:30, and descending-nighttime, ~01:30 overpasses), and that of Terra is on ~22:30 and ~10:30, respectively.We use MODIS Aqua and Terra average LST for GDD from MODIS.As a consequence, the satellite AGDD may be biased when compared with meteorological measurements, which in turn may limit the modeling of quantitative vegetation specific AGDD that drive phenophase transitions.
To characterize the seasonal progression of thermal time, we fitted the GDDs as a convex quadratic (CxQ) function of AGDD.The CxQ model have been successfully applied in temperate herbaceous vegetation and boreal ecosystems [21,47,65].
The intercept α is the start of observation period GDD value (which may not be zero due to the compositing), the linear parameter β affects the slope, and the quadratic parameter γ controls the curvature.Since our model is convex quadratic in shape, the sign of β is positive while the sign of γ is negative.
Two phenometrics were derived from the fitted parameter coefficients of the CxQ model.The peak height (PH) (Equation ( 4)) describes the maximum GDD from the fitted model, and the thermal time to peak (TTP) (Equation ( 5)) describes the amount of AGDD needed to reach the peak GDD: where α, β, and γ are the fitted parameter coefficients.
Our consideration of the growing season from day of year (DOY) 89-305 (roughly the beginning of April through October), allows us to characterize all the study cropland sites within a similar time frame.However, since our study sites have 12 • latitudinal range, there might be sites and years in the lower latitude winter cropland study sites that have above 0 • C temperature before DOY 89 or after DOY 305.Consequently, an underestimation of AGDD is possible for these sites or when either the early spring or the fall is warmer than normal.This warmer weather may affect the trajectory of NDVI or EVI as a function of AGDD, in particular for croplands with winter crops where the dormancy starts after DOY 305 and ends before DOY 89 of the following year, possibly biasing the TTP GDD .However, other factors such as different sowing dates, management practices, changes in crops under production, and field conditions around sowing and harvesting, may also affect cropland phenology.
We calculated the deviation of individual years' GDD time series from the CxQ model based on the average of eight years of GDD data.The root mean square deviation (RMSD) was used to calculate the average magnitude of the annual GDD deviation.Although eight years' time span is too a short period to construct a climatology, the expectation based on the shorter time series still provides a useful contrast to the highlight seasonal differences in particular years.
NDVI and EVI were calculated from the MODIS NBAR reflectance data using formulas developed by Tucker [8] and Huete et al. [18], respectively.We have plotted time series of NDVI and EVI (from the VNIR MODIS) as a function of AGDD (from the passive microwave AMSR-E).In the result section, we have described AGDDs at peak VIs based on growing season VI maximum values.AGDDs at 95% of the first peak VIs (peak height determined using growing season maximum value) were used to analyze temporal and spatial LSP responses to their thermal regime.In each AMSR-E pixel there are 25 to 36 corresponding MODIS pixels, depending on the latitude of the site.We generated corresponding MODIS and AMSR-E AGDDs for the 95% value of the initial seasonal peaks of NDVI and EVI.Since the peak VIs were selected based on the seasonal maximum value, sometimes these single peak values may occur in error.To avoid such errors, we use the 95% value of these peak VIs.While the NDVI, EVI and AGDD for the MODIS data were generated for each MODIS pixel within an AMSR-E pixel, the same AMSR-E AGDD was used for all the MODIS pixels within the same AMSR-E pixel.We assume that air temperature did not significantly vary within the AMSR-E pixel, an assumption that is reasonable in relatively flat terrain that is distant from large bodies of water.
Here we compare the seasonal development of MODIS VIs at the MODIS Climate Modeling Grid (CMG) resolution (0.05 • ) as a function of the AGDD derived from MODIS and from AMSR-E AGDDs.

Land Surface Seasonalities of Growing Degree-Days
In this subsection, we present the comparisons of satellite GDDs with meteorological station GDDs for quality assessment.The convex quadratic (CxQ) models of GDD are followed.We then present the AMSR-E thermal time to peak (TTP) (AGDD) in relation to latitude.Finally, AMSR-E GDD time series deviations from the multi-year average GDD model is presented.

Comparisons of GDD Derived from Satellite and Station Data
Meteorological stations 10 km buffer zone displayed a similar cropland (CRP) dominated average land use land cover for 2003-2010 study period (>90% cropland; Table 4).Saratov 4 is an AMSR-E pixel study site with considerable proportion of crop-natural vegetation mosaic (CNVM; 12%), while Kharkiv_2 meteorological station buffer has 7% urban and built-up (UBU) area.Comparison between Satellite GDD and station GDD was made using the 2003 GDD as nearly all the selected stations (that are within 25 km distance from the centers of AMSR-E pixel study sites) had full data during this year's growing period.Root mean square errors/deviations (RMSE or RMSD) between AMSR-E GDD and station GDD were lower than those between the MODIS GDD and station GDD (14-24 • C versus 16-30 • C; Table 5).Satellite GDD from both AMSR-E and MODIS showed strong linear correspondence with meteorological station GDD with r 2 > 0.87 (the covariate station GDD is correlated with response variable Satellite GDD; Figure 3 and Table 5).Overall, r 2 from the station and AMSR-E GDD fits was higher than those between the station and MODIS GDDs.GDDs from both satellite sensors had relatively higher magnitude compared to station GDDs.Based on analysis of covariance (ANCOVA), slopes from the linear regression fit of the two satellite sensors GDD with station GDD were partially significantly different (three out of six sites).That means that the covariate (station GDD) is not correlated to the independent variables (factor; sensors) for three sites while it is correlated for the rest three sites at α = 95% significance level.In addition, the ANCOVA model revealed that intercepts resulted from the linear regression fit between station GDD and satellite GDD are not significantly different by sensor in four sites (p > 0.05) while it is significantly different for two sites (p ≤ 0.05).Further analysis that incorporates many station data is needed to reach at a conclusion that linear fit of station GDD with AMSR-E GDD is significantly different from that with the MODIS GDD.We further closely examined the correspondence between daily station GDDs and AMSR-E GDDs.AMSR-E GDD not only tracked the seasonal dynamics of the station GDDs, but it also tracked the finer weekly and monthly GDD dynamics (Figure 4 left).The root mean squared difference (RMSD) between these daily datasets from the two sources was very low (1.86-3.16°C).This result is in line with [57] that found a bias of 1.0-3.5 K between the 2003 WMO weather station and AMSR-E daily air temperature data on vegetated lands in the Northern Hemisphere.It is also important to    We further closely examined the correspondence between daily station GDDs and AMSR-E GDDs.AMSR-E GDD not only tracked the seasonal dynamics of the station GDDs, but it also tracked the finer weekly and monthly GDD dynamics (Figure 4 left).The root mean squared difference (RMSD) between these daily datasets from the two sources was very low (1.86-3.16°C).This result is in line with [57] that found a bias of 1.0-3.5 K between the 2003 WMO weather station and AMSR-E daily air temperature data on vegetated lands in the Northern Hemisphere.It is also important to We further closely examined the correspondence between daily station GDDs and AMSR-E GDDs.AMSR-E GDD not only tracked the seasonal dynamics of the station GDDs, but it also tracked the finer weekly and monthly GDD dynamics (Figure 4 left).The root mean squared difference (RMSD) between these daily datasets from the two sources was very low (1.86-3.16• C).This result is in line with [57] that found a bias of 1.0-3.5 K between the 2003 WMO weather station and AMSR-E daily air temperature data on vegetated lands in the Northern Hemisphere.It is also important to note that there are times that the magnitude of GDDs from the two data sources is different, despite the similarities in dynamics, and this divergence occurs earlier and later in the growing season (e.g. Figure 4 left).These GDDs display strong linear correspondence with r 2 > 0.82 (e.g., Figure 4 right).However, there is a bias in these correspondence (as shown by the intercept) and underestimation (displayed by the slope) of the AMSR-E GDD relative to that of the station.

Convex Quadratic Models of GDD
The CxQ model fit well the relationship between GDD and AGDD for both the AMSR-E ta and MODIS LST, with coefficients of determination (r 2 ) ranging from 0.88 to 0.98.In contrast to the traditional approach of fitting GDD as a function of day of year (DOY), the CxQ fits for GDD as a function of AGDD displayed clear latitudinal shifts and more finely resolved intra-seasonal dynamics of GDD (Figure 5a,c versus Figure 5b,d).The more northern sites displayed a shorter summer growing season with narrower base and shorter peak GDD curves (Figure 5b), while the more southern sites displayed a longer growing season with broader base and larger peak GDD curves (Figure 5d).The fits for the two datasets (AMSR-E and MODIS) exhibited strong similar seasonalities (Figure 5).However, in 76% of our study sites, AMSR-E AGDD is larger than MODIS AGDD even though this difference is not statistically significant.The discrepancy may result from one or more of the following factors: (1) differences between surface air temperature and land surface temperature; (2) spatial scales of the retrievals-a single 25 km pixel versus ~25 pixels of 5.6 km; (3) differences between accumulating daily GDDs over eight days (AMSR-E) versus an eight-day composite (MODIS); and (4) residual atmospheric contamination in the MODIS retrievals.note that there are times that the magnitude of GDDs from the two data sources is different, despite the similarities in dynamics, and this divergence occurs earlier and later in the growing season (e.g. Figure 4 left).These GDDs display strong linear correspondence with r 2 > 0.82 (e.g., Figure 4 right).However, there is a bias in these correspondence (as shown by the intercept) and underestimation (displayed by the slope) of the AMSR-E GDD relative to that of the station.

Convex Quadratic Models of GDD
The CxQ model fit well the relationship between GDD and AGDD for both the AMSR-E ta and MODIS LST, with coefficients of determination (r 2 ) ranging from 0.88 to 0.98.In contrast to the traditional approach of fitting GDD as a function of day of year (DOY), the CxQ fits for GDD as a function of AGDD displayed clear latitudinal shifts and more finely resolved intra-seasonal dynamics of GDD (Figure 5a,c versus Figure 5b,d).The more northern sites displayed a shorter summer growing season with narrower base and shorter peak GDD curves (Figure 5b), while the more southern sites displayed a longer growing season with broader base and larger peak GDD curves (Figure 5d).The fits for the two datasets (AMSR-E and MODIS) exhibited strong similar seasonalities (Figure 5).However, in 76% of our study sites, AMSR-E AGDD is larger than MODIS AGDD even though this difference is not statistically significant.The discrepancy may result from one or more of the following factors: (1) differences between surface air temperature and land surface temperature; (2) spatial scales of the retrievals-a single 25 km pixel versus ~25 pixels of 5.6 km; (3) differences between accumulating daily GDDs over eight days (AMSR-E) versus an eight-day composite (MODIS); and (4) residual atmospheric contamination in the MODIS retrievals.Generally, higher and middle latitude study sites (>49 • N) have higher AMSR-E AGDD while those at lower latitudes have higher MODIS AGDD.This discrepancy may arise from the differences in LST versus surface air temperature, particularly during periods of low vegetation cover.All 49 study sites displayed consistently good fits with the CxQ model with coefficients of determination ranging from 0.88 to 0.98 (Table S1).The CxQ models fitted to the AMSR-E surface air temperature data displayed higher coefficients of determination for almost every site (42 of 49 sites) compared to the MODIS LST data (Table S1) and, considered as a group, the AMSR-E model fits were significantly (p < 0.00001) better than the MODIS model fits.
GDD at the beginning and end of the observational season were well above 0 • C at lower latitude sites (Figure 5c,d), but closer to zero at higher latitude sites (Figure 5a,b).The end-of-season AGDD was also larger (>4000 • C) at the lower latitude sites (Figure 5d) and smaller (~3000 • C) at the higher latitude sites (Figure 5b).For the rest of the analyses, we used the GDD and AGDD calculated from the AMSR-E surface air temperature data (unless stated otherwise).

Latitudinal Distributions of Thermal Regimes
Average daily GDD for 2003 through 2010 exhibited a strong inverse linear relationship with respect to latitude with r 2 > 0.90, in accordance with the long-established geographic fact that, at comparable elevation, temperature decreases as the absolute value of latitude increases (Figure 6a).The Thermal Time to Peak (TTP GDD ) graphed as a function of latitude revealed strong significant negative linear relationship (Figure 6b).The coefficient of determination for the sites in Ukraine was 0.81 and, for the sites in Russia and Kazakhstan combined, it is 0.92.Saratov 1 (RU) tends to be an outlier with higher TTP, likely due to its location farther inland from the Volga River.Across the 12 • latitudinal range of the study sites, there is a maximum difference in day length of two hours at the summer solstice.Generally, higher and middle latitude study sites (>49°N) have higher AMSR-E AGDD while those at lower latitudes have higher MODIS AGDD.This discrepancy may arise from the differences in LST versus surface air temperature, particularly during periods of low vegetation cover.All 49 study sites displayed consistently good fits with the CxQ model with coefficients of determination ranging from 0.88 to 0.98 (Table S1).The CxQ models fitted to the AMSR-E surface air temperature data displayed higher coefficients of determination for almost every site (42 of 49 sites) compared to the MODIS LST data (Table S1) and, considered as a group, the AMSR-E model fits were significantly (p < 0.00001) better than the MODIS model fits.
GDD at the beginning and end of the observational season were well above 0 °C at lower latitude sites (Figure 5c,d), but closer to zero at higher latitude sites (Figure 5a,b).The end-of-season AGDD was also larger (>4000 °C) at the lower latitude sites (Figure 5d) and smaller (~3000 °C) at the higher latitude sites (Figure 5b).For the rest of the analyses, we used the GDD and AGDD calculated from the AMSR-E surface air temperature data (unless stated otherwise).

Latitudinal Distributions of Thermal Regimes
Average daily GDD for 2003 through 2010 exhibited a strong inverse linear relationship with respect to latitude with r 2 > 0.90, in accordance with the long-established geographic fact that, at comparable elevation, temperature decreases as the absolute value of latitude increases (Figure 6a).The Thermal Time to Peak (TTPGDD) graphed as a function of latitude revealed strong significant negative linear relationship (Figure 6b).The coefficient of determination for the sites in Ukraine was 0.81 and, for the sites in Russia and Kazakhstan combined, it is 0.92.Saratov 1 (RU) tends to be an outlier with higher TTP, likely due to its location farther inland from the Volga River.Across the 12° latitudinal range of the study sites, there is a maximum difference in day length of two hours at the summer solstice.

Time Series of GDD Residuals
During warmer years, GDDs were well above the GDDs predicted by the multi-year average model, yielding positive GDD residuals and a larger root mean squared difference (RMSD) (Figure 7, right; RMSD = 46.3°C).The 2010 Russian heat wave that devastated grain production in European Russia [70], had a clear impact at Orenburg (Figure 7).During cooler years, GDDs were below the predicted GDDs yielding negative GDD residuals and a larger RMSD (Figure 7, left; RMSD = 22.3 °C).There were also some years that had GDDs that were comparable with the GDDs predicted by the average model, resulting in a smaller RMSD (Figure 7, center; RMSD = 19.5 °C).

Time Series of GDD Residuals
During warmer years, GDDs were well above the GDDs predicted by the multi-year average model, yielding positive GDD residuals and a larger root mean squared difference (RMSD) (Figure 7, right; RMSD = 46.3• C).The 2010 Russian heat wave that devastated grain production in European Russia [70], had a clear impact at Orenburg (Figure 7).During cooler years, GDDs were below the predicted GDDs yielding negative GDD residuals and a larger RMSD (Figure 7, left; RMSD = 22.3 • C).There were also some years that had GDDs that were comparable with the GDDs predicted by the average model, resulting in a smaller RMSD (Figure 7, center; RMSD = 19.5

Land Surface Phenologies of Vegetation Indices
In this subsection we first present the general LSP characteristics of the NDVI and EVI time series, then describe the cropland dynamics using the VIs in space and time, and we conclude with the responses of peak VIs to thermal time measured by AMSR-E and MODIS sensors.

NDVI and EVI
Synergistic use of the MODIS vegetation indices (VIs) NDVI and EVI with AMSR-E AGDD captures well the LSP typical for crops in the mid-to-high latitudes that are planted in spring and harvested in the fall (Figure 8).Cropland NDVI and EVI start to rise (green-up onset) as crops start to grow in spring, ascend to their peak seasonal value (maturity), and then decline during senescence and crop drying before harvest.Winter crops are planted in the fall and grow before a period of dormancy over the winter, then emerge in the spring and rapidly develop, before crop drying and harvesting in early summer.Both spring and winter crops could be found with the study region, but winter crops are restricted to the south tier of the region in Ukraine and in southern Russia north of the Black Sea.Since vegetation in the mid-latitudes is temperature-limited, the time series of each VI as a function of AGDD can be well fitted by the CxQ model [36,48].A general unimodal pattern of LSP was evident across the study region.However, there were variations to this pattern, as detailed below.

Land Surface Phenologies of Vegetation Indices
In this subsection we first present the general LSP characteristics of the NDVI and EVI time series, then describe the cropland dynamics using the VIs in space and time, and we conclude with the responses of peak VIs to thermal time measured by AMSR-E and MODIS sensors.

NDVI and EVI
Synergistic use of the MODIS vegetation indices (VIs) NDVI and EVI with AMSR-E AGDD captures well the LSP typical for crops in the mid-to-high latitudes that are planted in spring and harvested in the fall (Figure 8).Cropland NDVI and EVI start to rise (green-up onset) as crops start to grow in spring, ascend to their peak seasonal value (maturity), and then decline during senescence and crop drying before harvest.Winter crops are planted in the fall and grow before a period of dormancy over the winter, then emerge in the spring and rapidly develop, before crop drying and harvesting in early summer.Both spring and winter crops could be found with the study region, but winter crops are restricted to the south tier of the region in Ukraine and in southern Russia north of the Black Sea.Since vegetation in the mid-latitudes is temperature-limited, the time series of each VI as a function of AGDD can be well fitted by the CxQ model [36,48].A general unimodal pattern of LSP was evident across the study region.However, there were variations to this pattern, as detailed below.

Land Surface Phenologies of Vegetation Indices
In this subsection we first present the general LSP characteristics of the NDVI and EVI time series, then describe the cropland dynamics using the VIs in space and time, and we conclude with the responses of peak VIs to thermal time measured by AMSR-E and MODIS sensors.

NDVI and EVI
Synergistic use of the MODIS vegetation indices (VIs) NDVI and EVI with AMSR-E AGDD captures well the LSP typical for crops in the mid-to-high latitudes that are planted in spring and harvested in the fall (Figure 8).Cropland NDVI and EVI start to rise (green-up onset) as crops start to grow in spring, ascend to their peak seasonal value (maturity), and then decline during senescence and crop drying before harvest.Winter crops are planted in the fall and grow before a period of dormancy over the winter, then emerge in the spring and rapidly develop, before crop drying and harvesting in early summer.Both spring and winter crops could be found with the study region, but winter crops are restricted to the south tier of the region in Ukraine and in southern Russia north of the Black Sea.Since vegetation in the mid-latitudes is temperature-limited, the time series of each VI as a function of AGDD can be well fitted by the CxQ model [36,48].A general unimodal pattern of LSP was evident across the study region.However, there were variations to this pattern, as detailed below.

Types of Land Surface Phenologies in Croplands
The MODIS VIs as a function of AMSR-E AGDD revealed that lower latitude cropland sites (in southern Russia and Ukraine) had growing seasons that generally exhibited bimodal VI patterns.These sites are warmer and longer in terms of AGDD: first VI peaks at about AGDD = 1300 • C (Figure 9a) and AGDD = 900 • C (Figure 9c), and second VI peaks at about AGDD = 2800 • C (Figure 9a) and AGDD = 2400 • C (Figure 9c).In contrast, higher latitude cropland sites (northern sites in RU) had unimodal, shorter growing seasons in terms of AGDD (peak VIs at about AGDD = 1200 • C); Figure 9d).The middle latitude sites behave as unimodal but with a broader base of growing season (peak VI at about AGDD = 2000 • C; Figure 9b).Sites exhibiting bimodality in the VI curve are associated with the cultivation of both winter and spring crops in one year (Figure 9a,c).Those sites exhibiting unimodal VI curves were associated with either usually with only a spring crop (Figure 9b,d) or, less commonly, with only a winter crop (Figure 10, 2005-2009).

Types of Land Surface Phenologies in Croplands
The MODIS VIs as a function of AMSR-E AGDD revealed that lower latitude cropland sites (in southern Russia and Ukraine) had growing seasons that generally exhibited bimodal VI patterns.These sites are warmer and longer in terms of AGDD: first VI peaks at about AGDD = 1300 °C (Figure 9a) and AGDD = 900 °C (Figure 9c), and second VI peaks at about AGDD = 2800 °C (Figure 9a) and AGDD = 2400 °C (Figure 9c).In contrast, higher latitude cropland sites (northern sites in RU) had unimodal, shorter growing seasons in terms of AGDD (peak VIs at about AGDD = 1200 °C); Figure 9d).The middle latitude sites behave as unimodal but with a broader base of growing season (peak VI at about AGDD = 2000 °C; Figure 9b).Sites exhibiting bimodality in the VI curve are associated with the cultivation of both winter and spring crops in one year (Figure 9a,c).Those sites exhibiting unimodal VI curves were associated with either usually with only a spring crop (Figure 9b,d  Winter grains grown in the region-wheat, rye, and barley-are planted and emerge in the fall before becoming dormant over the winter.Resuming growth in the early spring when the air and soil temperatures warm, they produce a rapid increase in VI earlier in the year.Winter grains' ripening and senescence phase occurs in late spring through early summer, expressed by the rapid VI decrease by mid-summer.A second increase in VI is due to late spring planting and emergence with the cycle peaking late summer with early fall harvest.For areas that grow winter grains, another smaller VIs peak appeared in mid-fall (Figure 9a), corresponding to the emergence and growth of the next year's crop before winter dormancy.Winter grains grown in the region-wheat, rye, and barley-are planted and emerge in the fall before becoming dormant over the winter.Resuming growth in the early spring when the air and soil temperatures warm, they produce a rapid increase in VI earlier in the year.Winter grains' ripening and senescence phase occurs in late spring through early summer, expressed by the rapid VI decrease by mid-summer.A second increase in VI is due to late spring planting and emergence with the cycle peaking late summer with early fall harvest.For areas that grow winter grains, another smaller VIs peak appeared in mid-fall (Figure 9a), corresponding to the emergence and growth of the next year's crop before winter dormancy.

Interannual Variation in Land Surface Phenologies
Over the eight years of study, NDVI and EVI displayed similar patterns, but with clear seasonal and interannual variation.In the same location, some years displayed growing season bimodality while others exhibited unimodality.Variation in the timing of peak VI was higher in some locations than others.For example, Simferopol, located in Crimea, had three bimodal years (2003,2004 Some sites exhibited changes in seasonal pattern during the study period.We grouped the phenological patterns into three broad categories: no change, one change, or two changes (Figure 11 and Table S2).Supplementary Figures S1-S49 illustrate the interannual variation of LSPs at each site.Those sites exhibiting only a single seasonal pattern (either unimodal or bimodal) throughout the study period were categorized as no change and totaled 33 sites (these were consistently unimodal seasonal pattern throughout the study period).The seven sites that changed from unimodal to

Interannual Variation in Land Surface Phenologies
Over the eight years of study, NDVI and EVI displayed similar patterns, but with clear seasonal and interannual variation.In the same location, some years displayed growing season bimodality while others exhibited unimodality.Variation in the timing of peak VI was higher in some locations than others.For example, Simferopol, located in Crimea, had three bimodal years (2003,2004 Some sites exhibited changes in seasonal pattern during the study period.We grouped the phenological patterns into three broad categories: no change, one change, or two changes (Figure 11 and Table S2).Supplementary Figures S1-S49 illustrate the interannual variation of LSPs at each site.Those sites exhibiting only a single seasonal pattern (either unimodal or bimodal) throughout the study period were categorized as no change and totaled 33 sites (these were consistently unimodal seasonal pattern throughout the study period).The seven sites that changed from unimodal to bimodal or from bimodal to unimodal were classified as one change.Finally, nine sites changes from one seasonal pattern to another and then back again (unimodal→bimodal→unimodal or bimodal→unimodal→bimodal) were classified as two changes.A geographic pattern of change in seasonal pattern is evident in Figure 11: only sites at the lower latitudes displayed any changes in seasonal pattern.Details are presented in Table S2.
Remote Sens. 2016, 8, 1016 16 of 26 bimodal or from bimodal to unimodal were classified as one change.Finally, nine sites changes from one seasonal pattern to another and then back again (unimodal→bimodal→unimodal or bimodal→unimodal→bimodal) were classified as two changes.A geographic pattern of change in seasonal pattern is evident in Figure 11: only sites at the lower latitudes displayed any changes in seasonal pattern.Details are presented in Table S2.Croplands in Ukraine exhibited a wider range of seasonal patterns: some sites changed from unimodal to bimodal (three sites), some from bimodal to unimodal (three sites), and some others from bimodal to unimodal and then back again (three sites; Table S2).Croplands in southern parts of Russia similarly experienced a cycling of seasonal patterns: bimodal to unimodal (one site), bimodal →unimodal→bimodal (three sites) and unimodal→bimodal→unimodal (three sites).All sites in Kazakhstan retained the unimodal seasonal pattern throughout the study period.
3.2.4.Sensor-Specific Differences in Thermal Time to Peak VI MODIS pixels within each AMSR-E pixel exhibit generally similar but distinct spatial and temporal thermal time patterns.MODIS and AMSR-E AGDDs at 95% of the initial peak VI values display a strong positive linear relationship across space and time with coefficients of determination ranging from 0.60 to 0.99 (Figure 12).This relationship was generally strong for the lower latitude sites compared to that of the higher latitudes.Although it could be expected that the VIs for each MODIS pixel are more related to the same pixel's LST.However, the strong correlation of the two sensors' AGDDs showed that the much coarser spatial resolution AMSR-E surface air temperature (ta) data did an excellent job of capturing the progress of finer spatial resolution MODIS VIs (Figure 12).Croplands in Ukraine exhibited a wider range of seasonal patterns: some sites changed from unimodal to bimodal (three sites), some from bimodal to unimodal (three sites), and some others from bimodal to unimodal and then back again (three sites; Table S2).Croplands in southern parts of Russia similarly experienced a cycling of seasonal patterns: bimodal to unimodal (one site), bimodal→unimodal→bimodal (three sites) and unimodal→bimodal→unimodal (three sites).All sites in Kazakhstan retained the unimodal seasonal pattern throughout the study period.
3.2.4.Sensor-Specific Differences in Thermal Time to Peak VI MODIS pixels within each AMSR-E pixel exhibit generally similar but distinct spatial and temporal thermal time patterns.MODIS and AMSR-E AGDDs at 95% of the initial peak VI values display a strong positive linear relationship across space and time with coefficients of determination ranging from 0.60 to 0.99 (Figure 12).This relationship was generally strong for the lower latitude sites compared to that of the higher latitudes.Although it could be expected that the VIs for each MODIS pixel are more related to the same pixel's LST.However, the strong correlation of the two sensors' AGDDs showed that the much coarser spatial resolution AMSR-E surface air temperature (ta) data did an excellent job of capturing the progress of finer spatial resolution MODIS VIs (Figure 12).bimodal or from bimodal to unimodal were classified as one change.Finally, nine sites changes from one seasonal pattern to another and then back again (unimodal→bimodal→unimodal or bimodal→unimodal→bimodal) were classified as two changes.A geographic pattern of change in seasonal pattern is evident in Figure 11: only sites at the lower latitudes displayed any changes in seasonal pattern.Details are presented in Table S2.Croplands in Ukraine exhibited a wider range of seasonal patterns: some sites changed from unimodal to bimodal (three sites), some from bimodal to unimodal (three sites), and some others from bimodal to unimodal and then back again (three sites; Table S2).Croplands in southern parts of Russia similarly experienced a cycling of seasonal patterns: bimodal to unimodal (one site), bimodal →unimodal→bimodal (three sites) and unimodal→bimodal→unimodal (three sites).All sites in Kazakhstan retained the unimodal seasonal pattern throughout the study period.
3.2.4.Sensor-Specific Differences in Thermal Time to Peak VI MODIS pixels within each AMSR-E pixel exhibit generally similar but distinct spatial and temporal thermal time patterns.MODIS and AMSR-E AGDDs at 95% of the initial peak VI values display a strong positive linear relationship across space and time with coefficients of determination ranging from 0.60 to 0.99 (Figure 12).This relationship was generally strong for the lower latitude sites compared to that of the higher latitudes.Although it could be expected that the VIs for each MODIS pixel are more related to the same pixel's LST.However, the strong correlation of the two sensors' AGDDs showed that the much coarser spatial resolution AMSR-E surface air temperature (ta) data did an excellent job of capturing the progress of finer spatial resolution MODIS VIs (Figure 12).(c,f)) for 95% of the initial peak NDVI (a-c) and EVI (d-f).There is strong linear relationship between AGDDs from AMSR-E and MODIS at 95% of initial peak VIs for these sites in 2003 with r 2 ranging from 0.88 to >0.99.This strong positive relationship was consistent in space and time with r 2 ranging from 0.60 to >0.99 across all study sites and years.

Shifting Seasonal Patterns: Real or Illusory?
How do we interpret the apparent variation between double and single seasonal patterns that were evident in southern Ukraine and southern Russia?Note that the southern tier of the study region hosts winter crops and can also support spring crops.In Figure 10, fall green-up is evident in 2007 and 2008 due to winter crops.The 2003 European heat wave affected winter croplands in western and southern Ukraine [71], while the 2007 Ukrainian heat wave devastated spring grains in 2007 [72].These impacts are evident in the VI curves of Figure 10.Figures 8-10 give representative examples of land surface phenologies exhibiting unimodality and bimodality with variation in the placement of the single and double modes.Should we interpret a unimodal LSP as a single annual crop?Clearly, the answer is not always.The more northerly sites in Russia and Kazakhstan exhibit a single mode that peaks roughly mid-season (blue circles in Figure 11).These sites do have a single annual spring crop.In contrast, there are multiple instances of early season peaks, indicative of a winter crop, with no second peak, especially in the southern tier of the study region.These patterns-such as seen for 2005-2009 in Figure 10-could indicate either changes in cultivation practices or crop failures.
A study of cropping frequency in European Russia found areas in southern Russia, e.g., Stavropol (site 2 in this study) and Saratov (sites 21 and 25) were successfully cropped during fewer years than areas located farther north in the grain belt, e.g., Cheboksary and Kazan (sites 48 and 49) [73].Lower cropping frequency was a result of longer fallowing periods and hotter, drier summer conditions.However, they also reported widespread shifts from cereals to soybean, sunflower, and rapeseed as well as changes in crop rotation periods, including changes in fallow periods [73,74].Changes have been occurring across this vast agricultural landscape, but the proximate causes of those changes may not be readily accessible to remote sensing [20,[75][76][77][78][79][80].However, there is a rich set of variations in land surface phenologies across the study region through the study period (cf.Supplementary Figures S1-S49) that the AGDD based on the AMSR-E surface temperature data makes more readily apparent.

Growing Degree-Day Residuals versus Vegetation Indices
Rainfed croplands at mid-latitude cropland are thermally limited [81]; thus, vegetation dynamics are strongly linked to the seasonal progression of temperature.Residuals resulting from How do we interpret the apparent variation between double and single seasonal patterns that were evident in southern Ukraine and southern Russia?Note that the southern tier of the study region hosts winter crops and can also support spring crops.In Figure 10, fall green-up is evident in 2007 and 2008 due to winter crops.The 2003 European heat wave affected winter croplands in western and southern Ukraine [71], while the 2007 Ukrainian heat wave devastated spring grains in 2007 [72].These impacts are evident in the VI curves of Figure 10.Figures 8-10 give representative examples of land surface phenologies exhibiting unimodality and bimodality with variation in the placement of the single and double modes.Should we interpret a unimodal LSP as a single annual crop?Clearly, the answer is not always.The more northerly sites in Russia and Kazakhstan exhibit a single mode that peaks roughly mid-season (blue circles in Figure 11).These sites do have a single annual spring crop.In contrast, there are multiple instances of early season peaks, indicative of a winter crop, with no second peak, especially in the southern tier of the study region.These patterns-such as seen for 2005-2009 in Figure 10-could indicate either changes in cultivation practices or crop failures.
A study of cropping frequency in European Russia found areas in southern Russia, e.g., Stavropol (site 2 in this study) and Saratov (sites 21 and 25) were successfully cropped during fewer years than areas located farther north in the grain belt, e.g., Cheboksary and Kazan (sites 48 and 49) [73].Lower cropping frequency was a result of longer fallowing periods and hotter, drier summer conditions.However, they also reported widespread shifts from cereals to soybean, sunflower, and rapeseed as well as changes in crop rotation periods, including changes in fallow periods [73,74].Changes have been occurring across this vast agricultural landscape, but the proximate causes of those changes may not be readily accessible to remote sensing [20,[75][76][77][78][79][80].However, there is a rich set of variations in land surface phenologies across the study region through the study period (cf.Supplementary Figures S1-S49) that the AGDD based on the AMSR-E surface temperature data makes more readily apparent.

Growing Degree-Day Residuals versus Vegetation Indices
Rainfed croplands at mid-latitude cropland are thermally limited [81]; thus, vegetation dynamics are strongly linked to the seasonal progression of temperature.Residuals resulting from the comparison of average GDD as a convex quadratic function of AGDD with average GDD time series showed clear patterns associated with average VI time series for 2003-2010 across sites.Negative GDD residuals occurred at times of higher VIs, presumably due to higher latent heat flux; whereas, positive GDD residuals corresponded to the times earlier and later in the growing season when VIs were low, presumably due to higher sensible heat flux.
For example, Figure 13a shows the average NDVI and average GDD residuals along with maximum and minimum values as error bars for six cropland sites and two cropland/natural vegetation mosaic sites both in Kazakhstan at a similar latitude.These sites show clear unimodality with the cropland/natural vegetation mosaic exhibiting distinctly higher VI earlier and later than in the croplands.However, the GDD residual patterns are comparable but have distinctive responses to a given heat accumulation.Figure 13b,c display the average NDVI and GDD residuals for the southernmost two sites that experienced bimodal and unimodal cropland pattern alternately: Figure 13b,c shows the average for the bimodal (unimodal) years.The bimodal seasonal pattern results in a more complicated relationship between NDVI and the GDD residuals compared to the sites with unimodal seasonal pattern.dominated surface flux due to canopy development back again to sensible heat dominance as the canopy fades later season.Note that these are not direct observations of the surface but, rather, the difference between the model of average GDD as a quadratic function of AGDD and the average AMSR-E GDD observations.The information embedded in the model lack of fit arises from changes in the surface condition, viz., transpiring vegetation, which affects the surface air temperature.That such a pattern can be sensed at the coarse spatial resolution of AMSR-E grid cells is remarkable and points to the need for additional work to explore how much additional information can be gleaned from these time series.However, to stress again, these informative residuals result from "filtering" the time series using a CxQ model of GDD, which is, in turn, motivated by the widespread pattern observed in mid-latitude landscapes dominated by herbaceous vegetation (e.g., grasslands and croplands) that the surface air temperature is primarily driven by insolation [57].
We have shown that the time series of residuals between observed GDDs and the GDDs predicted from the CxQ model can show distinct shifts in air temperature corresponding to the growth and development and eventual senescence of the vegetated canopy.Specifically, above average air temperatures (positive residuals) occur in the weeks prior to increasing vegetation indices, followed by below average temperatures (negative residuals) as the vegetation indices are near peak, and then back to above average temperatures during vegetation senescence.This sequence points to the detection of shifts in the surface energy balance from high values of the Bowen ratio to low and back again to high.More study is needed to see how far this provocative finding can be pushed.Is it apparent over land covers other than croplands or in other locations?How does it compare to both coarser and finer remote sensing data streams?

Heat Wave Impacts on LSPs
The influence of two regional heat waves-the 2010 heat wave in Russia and Kazakhstan [70,82,83] and the 2007 heat wave in Ukraine [72]-were evident in LSP patterns generated from the MODIS and AMSR-E data.Western Russia experienced an extraordinary heat wave in 2010, in which July was the warmest month since 1980 and numerous locations setting all-time maximum temperature records [70].The 2003 European heat wave affected western and southern Ukraine but was only slightly evident in the seasonal GDD pattern and not in its accumulation [71].Winter grains of croplands in those areas appeared affected, displaying lower VI peaks for winter grains but higher VI peaks for spring grains.
For example, Figure 14 presented sample areas that combine GDD and VIs.GDD and VIs for the 2010 and 2007 heat wave years and average years for the same place are presented.During the These residuals appear to track the shift from sensible heat dominated surface flux to latent heat dominated surface flux due to canopy development back again to sensible heat dominance as the canopy fades later season.Note that these are not direct observations of the surface but, rather, the difference between the model of average GDD as a quadratic function of AGDD and the average AMSR-E GDD observations.The information embedded in the model lack of fit arises from changes in the surface condition, viz., transpiring vegetation, which affects the surface air temperature.That such a pattern can be sensed at the coarse spatial resolution of AMSR-E grid cells is remarkable and points to the need for additional work to explore how much additional information can be gleaned from these time series.However, to stress again, these informative residuals result from "filtering" the time series using a CxQ model of GDD, which is, in turn, motivated by the widespread pattern observed in mid-latitude landscapes dominated by herbaceous vegetation (e.g., grasslands and croplands) that the surface air temperature is primarily driven by insolation [57].
We have shown that the time series of residuals between observed GDDs and the GDDs predicted from the CxQ model can show distinct shifts in air temperature corresponding to the growth and development and eventual senescence of the vegetated canopy.Specifically, above average air temperatures (positive residuals) occur in the weeks prior to increasing vegetation indices, followed by below average temperatures (negative residuals) as the vegetation indices are near peak, and then back to above average temperatures during vegetation senescence.This sequence points to the detection of shifts in the surface energy balance from high values of the Bowen ratio to low and back again to high.More study is needed to see how far this provocative finding can be pushed.Is it apparent over land covers other than croplands or in other locations?How does it compare to both coarser and finer remote sensing data streams?

Heat Wave Impacts on LSPs
The influence of two regional heat waves-the 2010 heat wave in Russia and Kazakhstan [70,82,83] and the 2007 heat wave in Ukraine [72]-were evident in LSP patterns generated from the MODIS and AMSR-E data.Western Russia experienced an extraordinary heat wave in 2010, in which July was the warmest month since 1980 and numerous locations setting all-time maximum temperature records [70].The 2003 European heat wave affected western and southern Ukraine but was only slightly evident in the seasonal GDD pattern and not in its accumulation [71].Winter grains of croplands in those areas appeared affected, displaying lower VI peaks for winter grains but higher VI peaks for spring grains.
For example, Figure 14 presented sample areas that combine GDD and VIs.GDD and VIs for the 2010 and 2007 heat wave years and average years for the same place are presented.During the average years, peak NDVIs and EVIs were above 0.60 and 0.35, respectively, while these peaks were attained at about 1230 • C AGDD for Kuybuskev 2 (Russia) and 960 • C AGDD for Mykolayive (Ukraine) (Table 6).Peak GDD for these sites was below 200 • C, which is equivalent to a daily average of ~25 • C over 8 days).In contrast, during the heat wave years, the maximum attained NDVIs and EVIs were below 0.50 and 0.25, respectively, while these maxima were attained earlier at about 960 • C AGDD for Kuybuskev 2 and 660 • C AGDD for Mykolayive.Peak GDD for these sites was at about 250 • C (≈31 • C average daily GDD), well above the average years' peak GDD.Crops were devastated by the heat waves before maturity resulting in crop production very much reduced compared to cooler, wetter years [84].
We grouped sites by the effects of the 2007 and 2010 major heat waves into four categories: (a) affected only in 2010; (b) affected only in 2007; (c) affected in both years; and (d) affected in neither year (Figure 15).The northern portion of the study region (both Russia and Kazakhstan) was strongly affected by the exceptional 2010 heat wave, southern Russia and eastern Ukraine croplands were affected by both heat waves, the southern portion of Ukraine was affected only by the 2007 heat wave, and some parts of western Ukraine, eastern and southern Russia were affected by neither heat wave.The heat wave affected areas exhibited AGDD time series were consistently higher than the other years (data not shown).Crop production was devastated in areas that were affected by these heat waves [72,82].

Conclusions
Although the spatial resolution of AMSR-E data appears coarse relative to imaging sensors, having a time series of surface air temperature every 25 km translates into a significant increase in spatial coverage over the weather station networks in most countries in the world, including developing nations that do not have dense networks of meteorological stations.Moreover, due to the long wavelengths of the microwaves, cloudiness does not pose the problem it does at the much shorter wavelengths of the thermal infrared.Instead, data gaps arise over frozen surfaces, during very heavy precipitation events, as a result of radio frequency interference over some parts of the planet, and due to orbital paths.Simple retrospective smoothing of the AMRS-E data yielded air temperature time series that compared well with ground observations.In the past AGDD has been calculated from meteorological station data (point coverage), meteorological station reanalysis (moderate to coarse spatial resolution), model reanalysis (very coarse spatial resolution), and thermal remote sensing data (moderate spatial resolution but temporally composited).The balance of coarse spatial resolution with high temporal resolution, due to the availability of nighttime acquisitions and cloud penetration, makes the AMSR-E data attractive as an additional source of surface air temperature data.
We have demonstrated that the surface air temperature data retrieved from a passive microwave radiometer (AMSR-E) agrees well with the weekly, monthly, and seasonal oscillations of meteorological station air temperature data.AMSR-E GDD also showed lower RMSD from station GDD compared to the corresponding MODIS GDD.Both GDDs from AMSR-E and MODIS showed strong correlation with station air temperature GDD with slight superiority of the AMSR-E GDD, but with significantly different fit in some sites and not in some other.Further study using more station data is needed to conclude whether the agreement between the AMSR-E GDDs and the station GDDs significantly differs from agreement between the MODIS GDDs and the station GDDs.
We have demonstrated the temporal relationship between GDD and AGDD from AMSR-E could be parsimoniously represented using a simple parametric convex quadratic (CxQ) model, which is particularly helpful to portray a climatology.Moreover, this relationship has a clear geographic pattern in mid-latitude croplands: decreasing thermal time to peak GDD (TTP GDD ) as latitude increases.This concise parametric representation of the expected progression of thermal time may offer some advantages in modeling vegetation growth and development.
We have shown how the AGDD can be combined with vegetation indices from MODIS to illustrate land surface phenology (LSP).What is new here is not that the CxQ model works for LSP-that has been shown many times before [20,21,31,36,42,48,49,67]-but rather that the surface air temperature data from a passive microwave source can be used to calculate the AGDD, which can, in turn, be used to model VIs from the moderate spatial resolution sensor MODIS.
AMSR-E failed in October 2011.AMSR-2 is onboard the Global Change Observations Mission 1st-Water (GCOM-W1), renamed SHIZUKU.It was launched in May 2012 by the Japan Aerospace Exploration Agency (JAXA).Data products were released to the public in May 2013 (https://gcom-w1.jaxa.jp/).AMSR-2 is the improved version of AMSR-E that has similar functionality and data products.The NTSG group has continued the production of the enhanced land parameters using the AMSR-2 data streams (ftp://ftp.ntsg.umt.edu/pub/data/AMSR_Results/AMSR_E_2/) as well as filled the data gap using a Chinese microwave radiometer [85,86].Thus, the continued relevance of this approach is assured for the next several years, at least.
The failure of the radar in NASA's SMAP mission during the summer of 2015 has refocused attention on the L-band passive radiometer data that has a 36 km resolution.While L-band seems too long a wavelength to detect crop dynamics, it should work well for soil moisture monitoring as well as potentially providing insight on land surface dynamics in ecosystems dominated by woody species.The approach developed here may also find traction with the L-band data, if it can be blended with the AMSR-2 data.
Finally, although this paper has focused on the use of the air temperature data from AMSR-E, there is another variable from the NTSG's dataset that is relevant to vegetation dynamics, the vegetation

Figure 2 .
Figure 2. Study region cropland stability map superimposed with the 49 specific AMSR-E pixels selected for this study.The AMSR-E pixels are numbered by latitude starting from the most southern site.Name for each site is their closest large settlement (cf.Table3).Red squares are in Ukraine, cyan squares in Russia, and blue squares in Kazakhstan.

Figure 2 .
Figure 2. Study region cropland stability map superimposed with the 49 specific AMSR-E pixels selected for this study.The AMSR-E pixels are numbered by latitude starting from the most southern site.Name for each site is their closest large settlement (cf.Table3).Red squares are in Ukraine, cyan squares in Russia, and blue squares in Kazakhstan.

Figure 3 .
Figure 3. Scatter plots and linear regression fits of station GDD with satellite GDD-AMSR-E GDD (black circles) and MODIS GDD (blue diamonds)-at Simferopol, Ukraine (site 4) for 2003.The linear regression fit for the two datasets were high with r 2 of 0.93 for the AMSR-E GDD and r 2 of 0.89 for the MODIS GDD.

Figure 4 .
Figure 4. (left) Time series plot of AMSR-E (black) and station (red) daily GDD as a function of station AGDD for 2003 at Kirovohrad, Ukraine (site 15); and (right) linear regression fit of AMSR-E GDD with station GDD for the same dataset yielding strong correspondence with r 2 = 0.95.Note also the bias (intercept) and the underestimation (slope < 1) of the AMSR-E GDD relative to the station GDD.

Figure 3 .
Figure 3. Scatter plots and linear regression fits of station GDD with satellite GDD-AMSR-E GDD (black circles) and MODIS GDD (blue diamonds)-at Simferopol, Ukraine (site 4) for 2003.The linear regression fit for the two datasets were high with r 2 of 0.93 for the AMSR-E GDD and r 2 of 0.89 for the MODIS GDD.

Figure 3 .
Figure 3. Scatter plots and linear regression fits of station GDD with satellite GDD-AMSR-E GDD (black circles) and MODIS GDD (blue diamonds)-at Simferopol, Ukraine (site 4) for 2003.The linear regression fit for the two datasets were high with r 2 of 0.93 for the AMSR-E GDD and r 2 of 0.89 for the MODIS GDD.

Figure 4 .
Figure 4. (left) Time series plot of AMSR-E (black) and station (red) daily GDD as a function of station AGDD for 2003 at Kirovohrad, Ukraine (site 15); and (right) linear regression fit of AMSR-E GDD with station GDD for the same dataset yielding strong correspondence with r 2 = 0.95.Note also the bias (intercept) and the underestimation (slope < 1) of the AMSR-E GDD relative to the station GDD.

Figure 4 .
Figure 4. (left) Time series plot of AMSR-E (black) and station (red) daily GDD as a function of station AGDD for 2003 at Kirovohrad, Ukraine (site 15); and (right) linear regression fit of AMSR-E GDD with station GDD for the same dataset yielding strong correspondence with r 2 = 0.95.Note also the bias (intercept) and the underestimation (slope < 1) of the AMSR-E GDD relative to the station GDD.

Figure 6 .
Figure 6.Thermal climates as a function of latitude revealed by: (a) average daily GDD; and (b) Thermal Time to Peak (TTPGDD).Latitudes were the geographic centers of AMSR-E pixels.All 49 study sites are displayed in both figures.In (b), hollow circles = Russia, orange diamonds = Kazakhstan, and cyan crosses on blue background = Ukraine.Both panels show a general decrease in (a) GDD or (b) TTP as latitude increases from 44° to 56°N.The uppermost two hollow circles are the northernmost study sites (site 49, Kazan', Russia and site 48, Cheboksary, Russia), while the lowermost two hollow circles are the southernmost study sites (site 1, Cherkessk, Russia and 2, Stavropol, Russia).

Figure 6 .
Figure 6.Thermal climates as a function of latitude revealed by: (a) average daily GDD; and (b) Thermal Time to Peak (TTP GDD ).Latitudes were the geographic centers of AMSR-E pixels.All 49 study sites are displayed in both figures.In (b), hollow circles = Russia, orange diamonds = Kazakhstan, and cyan crosses on blue background = Ukraine.Both panels show a general decrease in (a) GDD or (b) TTP as latitude increases from 44 • to 56 • N. The uppermost two hollow circles are the northernmost study sites (site 49, Kazan', Russia and site 48, Cheboksary, Russia), while the lowermost two hollow circles are the southernmost study sites (site 1, Cherkessk, Russia and 2, Stavropol, Russia).

Figure 8 .
Figure 8.Average MODIS NDVI and EVI as a function of AMSR-E AGDD for Petropavlovsk 3, Kazakhstan (site 42) for 2003-2010.Note that the NDVI displays a larger dynamic range than the EVI.

Figure 8 .
Figure 8.Average MODIS NDVI and EVI as a function of AMSR-E AGDD for Petropavlovsk 3, Kazakhstan (site 42) for 2003-2010.Note that the NDVI displays a larger dynamic range than the EVI.

Figure 8 .
Figure 8.Average MODIS NDVI and EVI as a function of AMSR-E AGDD for Petropavlovsk 3, Kazakhstan (site 42) for 2003-2010.Note that the NDVI displays a larger dynamic range than the EVI.

Figure 10 .
Figure 10.NDVI and EVI interannual variability in one of the southernmost study sites (Simferopol, Ukraine at 45.6°N) from 2003-2010.Whether due to changes in cultivation practice or crop failures, the VI curves change from bimodal (2003-2004) to unimodal (2005-2009) and back to bimodal (2010).
, and 2010) and five unimodal years (2005-2009) (Figure 10).Variation was evident within bimodal growing seasons: in 2003 the second peak was higher (peak NDVI of 0.57 versus 0.48); in 2004 peaks were wider apart in time (difference between first and second peak AGDD of 1900 °C in 2004 versus 1550 °C in 2003); and in 2010 the second peak was much smaller (peak NDVI of 0.51 versus 0.59).Years with a unimodal peak generally had peaks that were larger compared to the bimodal peaks (peak NDVI of >0.63 versus <0.63).In 2007 and 2008 the peaks were earlier; while in 2003, 2009, and 2010 the peak heights were lower.Moreover, since Simferopol is a southern site, the growing season can start before DOY 89, but with interannual variation due to weather.The effect on the VI curves in Figure 10 is clear: the early season green-up pattern in 2003-2005 has initial NDVI (EVI) at or below 0.4 (0.2).In contrast, initial VIs in 2006-2010 are well above 0.4 and 0.2 for the NDVI and EVI, respectively.

Figure 10 .
Figure 10.NDVI and EVI interannual variability in one of the southernmost study sites (Simferopol, Ukraine at 45.6 • N) from 2003-2010.Whether due to changes in cultivation practice or crop failures, the VI curves change from bimodal (2003-2004) to unimodal (2005-2009) and back to bimodal (2010).
, and 2010) and five unimodal years (2005-2009) (Figure 10).Variation was evident within bimodal growing seasons: in 2003 the second peak was higher (peak NDVI of 0.57 versus 0.48); in 2004 peaks were wider apart in time (difference between first and second peak AGDD of 1900 • C in 2004 versus 1550 • C in 2003); and in 2010 the second peak was much smaller (peak NDVI of 0.51 versus 0.59).Years with a unimodal peak generally had peaks that were larger compared to the bimodal peaks (peak NDVI of >0.63 versus <0.63).In 2007 and 2008 the peaks were earlier; while in 2003, 2009, and 2010 the peak heights were lower.Moreover, since Simferopol is a southern site, the growing season can start before DOY 89, but with interannual variation due to weather.The effect on the VI curves in Figure 10 is clear: the early season green-up pattern in 2003-2005 has initial NDVI (EVI) at or below 0.4 (0.2).In contrast, initial VIs in 2006-2010 are well above 0.4 and 0.2 for the NDVI and EVI, respectively.

Figure 11 .
Figure 11.Changes in seasonal patterns by sites during 2003-2010: no change (blue circles), one change (white stars with pink borders), or two changes (red squares).Northern study sites displayed no change in seasonal pattern, while southern study sites experienced multiple changes.

Figure 11 .
Figure 11.Changes in seasonal patterns by sites during 2003-2010: no change (blue circles), one change (white stars with pink borders), or two changes (red squares).Northern study sites displayed no change in seasonal pattern, while southern study sites experienced multiple changes.

Figure 11 .
Figure 11.Changes in seasonal patterns by sites during 2003-2010: no change (blue circles), one change (white stars with pink borders), or two changes (red squares).Northern study sites displayed no change in seasonal pattern, while southern study sites experienced multiple changes.

Figure 12 .
Figure 12. (a-f) Comparison of AMSR-E AGDDta at 95% of the initial peak VIs with the MODIS AGDDlst at 95% of the initial peak VIs in 2003.The panels span a latitudinal gradient: lowest latitude (Cherkessk, Russia (a,d)), middle latitude (Sumy, Ukraine (b,e)), and highest latitude (Kazan', Russia(c,f)) for 95% of the initial peak NDVI (a-c) and EVI (d-f).There is strong linear relationship between AGDDs from AMSR-E and MODIS at 95% of initial peak VIs for these sites in 2003 with r 2 ranging from 0.88 to >0.99.This strong positive relationship was consistent in space and time with r 2 ranging from 0.60 to >0.99 across all study sites and years.

Figure 12 . 4 . 1 .
Figure 12. (a-f) Comparison of AMSR-E AGDDta at 95% of the initial peak VIs with the MODIS AGDDlst at 95% of the initial peak VIs in 2003.The panels span a latitudinal gradient: lowest latitude (Cherkessk, Russia (a,d)), middle latitude (Sumy, Ukraine (b,e)), and highest latitude (Kazan',Russia (c,f)) for 95% of the initial peak NDVI (a-c) and EVI (d-f).There is strong linear relationship between AGDDs from AMSR-E and MODIS at 95% of initial peak VIs for these sites in 2003 with r 2 ranging from 0.88 to >0.99.This strong positive relationship was consistent in space and time with r 2 ranging from 0.60 to >0.99 across all study sites and years.

Figure 13 .
Figure 13.Average GDD residuals and NDVI with error bars displaying maxima and minima for sites at similar latitude.(a) Sites 35-37 and 39-41 in Kazakhstan (NDVI-blue triangles and GDD residuals-black diamonds) and sites 34 and 42 in Kazakhstan (NDVI-green crosses and GDD residuals-gray circles); and (b,c) sites 1 and 2 in Russia (NDVI-blue diamonds and GDD residuals-black circles) for selected bimodal cropland pattern years and unimodal years, respectively.

Figure 15 .
Figure 15.Sites affected by: the 2010 heat wave (blue circles), the 2007 heat wave (red squares), both the 2010 and 2007 heat waves (white stars with pink borders), and sites affected by neither the 2010 nor the 2007 heat wave (cyan triangles).

Figure 15 .
Figure 15.Sites affected by: the 2010 heat wave (blue circles), the 2007 heat wave (red squares), both the 2010 and 2007 heat waves (white stars with pink borders), and sites affected by neither the 2010 nor the 2007 heat wave (cyan triangles).

Table 1 .
Percentage of days with meteorological station data missing during the growing season.

Table 2 .
Interpretative legend for Figure 1 that display IGBP MODIS 0.05 • land cover variation from 2003-2010 in the study region.The table shows how the color in the LC map (Figure

Table 3 )
. Red squares are in Ukraine, cyan squares in Russia, and blue squares in Kazakhstan.

Table 4 .
Average (2003Average ( -2010) MODIS land cover with the 10 km buffer zone for meteorological stations and corresponding AMSR-E pixel study sites.Bold indicates the dominant land cover (all sites are in crop dominated areas).

Table 5 .
Root mean square error (RMSE), linear regression model (intercept, slope and r 2 ), and analysis of covariance (ANCOVA) between station and satellite GDDs from AMSR-E and MODIS sensors.AMS = AMSR-E; MOD = MODIS.
* RMSE of observed satellite and meteorological station GDDs (not using fitted GDDs).

Table 5 .
Root mean square error (RMSE), linear regression model (intercept, slope and r 2 ), and analysis of covariance (ANCOVA) between station and satellite GDDs from AMSR-E and MODIS sensors.AMS = AMSR-E; MOD = MODIS.
* RMSE of observed satellite and meteorological station GDDs (not using fitted GDDs).

Table 5 .
Root mean square error (RMSE), linear regression model (intercept, slope and r 2 ), and analysis of covariance (ANCOVA) between station and satellite GDDs from AMSR-E and MODIS sensors.AMS = AMSR-E; MOD = MODIS.

Value) AMS MOD AMS MOD AMS MOD AMS MOD Slope
* RMSE of observed satellite and meteorological station GDDs (not using fitted GDDs).