Assessing the Impact of Climate Variability on Cropland Productivity in the Canadian Prairies Using Time Series MODIS FAPAR

Abstract

Cropland productivity is impacted by climate. Knowledge on spatial-temporal patterns of the impacts at the regional scale is extremely important for improving crop management under limiting climatic factors. The aim of this study was to investigate the effects of climate variability on cropland productivity in the Canadian Prairies between 2000 and 2013 based on time series of MODIS (Moderate Resolution Imaging Spectroradiometer) FAPAR (Fraction of Absorbed Photosynthetically Active Radiation) product. Key phenological metrics, including the start (SOS) and end of growing season (EOS), and the cumulative FAPAR (CFAPAR) during the growing season (between SOS and EOS), were extracted and calculated from the FAPAR time series with the Parametric Double Hyperbolic Tangent (PDHT) method. The Mann-Kendall test was employed to assess the trends of cropland productivity and climatic variables, and partial correlation analysis was conducted to explore the potential links between climate variability and cropland productivity. An assessment using crop yield statistical data showed that CFAPAR can be taken as a surrogate of cropland productivity in the Canadian Prairies. Cropland productivity showed an increasing trend in most areas of Canadian Prairies, in general, during the period from 2000 to 2013. Interannual variability in cropland productivity on the Canadian Prairies was influenced positively by rainfall variation and negatively by mean air temperature.

1. Introduction

Canada is one of the largest food producers and exporters in the global market (www.grainscanada.gc.ca). In addition to cultivar turnover and changes in agronomic management, variations of crop growth and productivity are mainly driven by climate variability [1,2,3,4]. Studies indicate that the climate is changing in Canada [5,6,7,8,9], and that it has a profound impact on crop production in the agriculture regions [2,3,5,10,11], especially in the Canadian Prairies that account for around 85% of Canada’s cropland area (http://www.statcan.gc.ca/). A better understanding of the impact of climate variables on cropland productivity can help improve crop management practices and the accuracy of crop production prediction, which is urgently needed with the increasing worldwide demand for food [12,13].

In the past few decades, two main approaches—statistical approach and crop model simulation—have been employed to investigate the influence of climatic variability on crop growth and productivity on the Canadian Prairies [3,8,10,14,15,16,17]. Based on long-term climate datasets, Qian et al. [5,11] analyzed the long-term variation of agro-climatic indices which are important for agricultural planning and management. They observed significant trends in agro-climatic indices (e.g., earlier start and delayed end of growing season) in Canada. Of all the agro-climatic factors, water-related factors or drought indices have played the most important role in variability in crop growth and productivity in the Canadian Prairies [2,3,8,18,19]. For instance, Mkhabela et al. [19] found that drought indices, especially water demand (evapotranspiration, ET) and the water balance index, can explain between 27% and 74% of the variation of wheat yield (p < 0.05) in the Canadian Prairies over the period 2003–2006. In addition, efforts to explore crop yield variation under future climate scenarios have been made using crop models and climate models [16,20,21]. Although these studies have greatly improved our understanding of the impact of climate change on crop growth, most of them were conducted based on data from sparse meteorological stations and annual surveyed crop yield. The results were unable to demonstrate the spatial heterogeneity associated with the relationship between climate indices and cropland productivity [22].

Vegetation indices and biophysical variables derived from remote sensing data, such as the Normalized Difference Vegetation index (NDVI) and the Fraction of Absorbed Photosynthetically Active Radiation (FAPAR), show a strong capability to monitor spatial-temporal variation of crop growth status and assess cropland productivity at the regional or global scale [4,12,13,23]. Over the past decades, a few long-term earth observation datasets have been made available, such as the TM/ETM+/OLI from the Landsat series satellites, the Advanced Very High Resolution Radiometer (AVHRR) from the NOAA satellites, the Moderate Resolution Imaging Spectroradiometer (MODIS) from the Terra/Aqua satellites, and the VEGETATION from the SPOT satellites [24,25]. These long-term datasets are useful in exploring how cropland productivity responds to climate variation (e.g., rainfall and temperature) and management practices [23,26]. For instance, Meroni et al. [4] found that crop biomass showed a strong positive correlation with growing season length (GSL) and a negative correlation with delays at the start of season (SOS) in Sahel. Based on time series NDVI, Tottrup and Rasmussen [26] suggested that the trend of crop yield in Senegal was not only affected by rainfall but was also significantly influenced by changes in land cover and land use. Milesi et al. [22] observed a declining trend in the relative growth rate of grain production in water-limited tropic countries (e.g., India) using the annual sum of the 12 monthly values of maximum NDVI. These studies indicated that long-term satellite data with their good spatial coverage could be useful in investigating the relationship between crop growth and climate variability. To date, in Canada, most studies have focused on estimating crop yield using remote sensing technology [12,13,27,28,29], while few studies have paid attention to the crop yield variation analysis using time series of remote sensing observations.

The main objectives of this study were: (1) to analyze the spatial and temporal variation of cropland productivity on the Canadian Prairies from 2000 to 2013 using satellite data; and (2) to explore the effects of growing season climate variables including temperature, precipitation, and radiation on the cropland productivity derived from satellite data.

2. Materials and Methods

2.1. FAPAR Product

FAPAR is a biophysical variable indicative of vegetation photosynthetic status [30,31,32,33]. It is directly related to primary productivity and is widely used in vegetation productivity modeling based on the light use efficiency theory [30,31]. Different from vegetation indices (VI), FAPAR has a clearer physical meaning [34], although some VIs (e.g., the enhanced vegetation index and the wide dynamic range vegetation index) have been found to have strong correlation with crop instantaneous productivity (e.g., GPP and NPP) [35,36]. Several studies have also demonstrated that cropland productivity is related more to FAPAR than to NDVI or leaf area index (LAI) [32,37,38]. FAPAR is also a critical component of the land-surface radiation budget that is useful in investigating the impact of climate variability on vegetation dynamics and has been recognized as one of the Essential Climate Variables (ECV) in the Global Terrestrial Observing System (GTOS) [39]. In this study, time series FAPAR from the MODIS FAPAR/LAI product (MOD15A2, collection 5) was employed to extract phenological metrics and derive cropland productivity. MODIS FAPAR is derived using physical models as the main algorithm and empirical relationships between NDVI and FAPAR as a backup algorithm when the physical models fail [40]. The product has a spatial resolution of 1 km with 8 days’ temporal composition to reduce cloud contamination. Validation studies have shown that it is sufficient for agricultural monitoring and cropland productivity estimation [41,42]. A total of 14 years’ MOD15A2 data from 2000 to 2013 were downloaded from the MODIS website (http://modis.gsfc.nasa.gov/). The product was processed using the MODIS re-projection tool (MRT, https://lpdaac.usgs.gov/tools/modis_reprojection_tool) to mosaic and extract FAPAR data.

2.2. Cropland Mask

To derive a general cropland mask for the period of 2000–2013, MODIS Collection 5 annual land cover product (MCD12Q1) from 2001 to 2012 was obtained from the MODIS website (http://modis.gsfc.nasa.gov/). The classification accuracy of the cropland type around the world was greater than 70% [43], which is satisfactory for the analyses of this study. In order to create an annual crop mask, two land cover types were extracted from the land cover type 1 based on the IGBP global vegetation classification scheme; cropland (code = 12) and cropland/natural vegetation mosaic (code = 14). The general cropland mask for the period of 2000–2013 was created by using any pixel that belonged to either code 12 or 14 in all years. As the resolution of MCD12Q1 is 500 m, the general cropland mask was resampled using the NEAREST method of the ArcGIS software to the same resolution as MODIS FAPAR. The cropland mask for Alberta (AB), Saskatchewan (SK) and Manitoba (MB) is shown in Figure 1a. The Canadian Prairie Provinces consists of the largest crop area in Canada and accounts for about 85% of the national total harvested area (http://www.statcan.gc.ca/). In general, the croplands on the Canadian Prairies are predominantly rainfed [44]. Canola, spring wheat, barley, and oats are the four major annual crops in this region [44]. The region is covered by three agro-climate zones; the sub-humid zone, the semi-arid zone and the arid zone (Figure 1b) [19]. The soils are dominated by Brown Chernozems in the arid zone, Dark-Brown Chernozems in the semi-arid zone and Black and Dark Gray Chernozems and Gray Luvisols in the sub-humid zone [13,45,46].

2.3. Climate Data

To assess the effects of climate factors on cropland productivity variation on the Canadian Prairies, two climate datasets were used: a daily 10-km gridded climate dataset generated by Natural Resources Canada (http://cfs.nrcan.gc.ca/projects/3/1), and the Modern Era Retrospective analysis for Research and Applications (MERRA) dataset from NASA (http://disc.sci.gsfc.nasa.gov/mdisc/). The 10-km gridded dataset contains daily maximum air temperature (°C), minimum air temperature (°C), and precipitation (mm) for the Canadian landmass south of 60°N between 2000 and 2013. The 10-km gridded data was previously interpolated from daily Environment Canada climate station observations using a thin plate smoothing spline surface fitting method implemented by ANUsplin V4.3 [17]. Daily mean temperature was calculated as the average of daily maximum and minimum temperatures. The MERRA dataset is a NASA 30-year (1979–present) reanalysis product using the Goddard Earth Observing System Data Assimilation System, Version 5 (GEOS-5), from which surface incident shortwave flux was obtained. Validations show good agreement between various meteorological factors (e.g., temperature, radiation, humidity, and energy balance) in the MERRA reanalysis dataset and surface meteorological data at the global scale [47,48]. The spatial resolution is 0.5° latitude by 0.67° longitude, and temporal resolution is three-hours.

2.4. Phenological Metrics

The Parametric Double Hyperbolic Tangent (PDHT) method, a model-fitting approach [4], was employed to derive the start (SOS) and the end of growing season (EOS) in this study. This approach is based on fitting the temporal profile of FAPAR with a mathematical model, and is useful in reducing the high-frequency noise induced by interference factors. The fitted time series of FAPAR was then used to derive SOS and EOS pixel by pixel on the cropland. SOS was defined as the day of year (DOY) at the local maxima of the second derivative from the rising side of the fitted curve, and EOS was defined as the DOY at the local maxima of the second derivative from the declining side of the fitted curve [49]. Cumulative FAPAR (CFAPAR) was calculated as the integration of the fitted daily FAPAR between SOS and EOS. To analyze the effect of climate variables on cropland productivity variation, the CFAPAR was resampled to 10-km spatial sampling interval to match that of the climate variables, using the NEAREST method of the ArcGIS software.

2.5. Crop Yield Statistics

Annual crop yield, harvest area and production data for the Canadian Prairies from 2000 to 2013 were obtained from Statistics Canada's socioeconomic database CANSIM (http://www5.statcan.gc.ca/cansim/). The crop data were collected through an annual farm sample survey program and are reported at the Census Agricultural Region (CAR) level for each crop type. CARs without complete yield records in each year were excluded from this study. The yield of barely, wheat, canola and oats in each CAR were then converted into net primary production (NPP, gC·m⁻²·yr⁻¹) based on the method provided by Prince, et al. [50] as follows,

$$NPP=\left(\left(Yield\times f_{mass}\right)/HI\times \left(1+RS\right)\right)\times f_{dry}\times f_{carbon}$$

(1)

where Yield is the crop yield obtained from the statistical data, HI is the harvest index, RS is the root to shoot ratio used to estimate the total biomass of the crop, f_mass is a factor to convert the yield in reporting units to a standard unit, f_dry is a factor to convert the mass to dry biomass, and f_carbon is a factor to convert the dry biomass to carbon (450 g·kg⁻¹). The conversion factors in different crop types were referenced from the study of Prince, et al. [50].

The final NPP integrated with different crop in each CAR was calculated on a unit per area basis (NPP_CAR, g·C·m⁻²) with the following equation:

$$NP{P}_{CAR}={{\displaystyle \sum }}_{i=0}^{n}NPP\left(i\right)\times Area\left(i\right)/{{\displaystyle \sum }}_{i=0}^{n}Area\left(i\right)$$

(2)

At the same time, the derived CFAPAR was averaged into FAPAR at the CAR level (FAPAR_CAR) based on the area of the cropland mask of the CAR. The correlation analysis between NPP_CAR and FAPAR_CAR was then conducted to evaluate the representativeness of CFAPAR for cropland productivity.

2.6. Statistical Analysis

2.6.1. Mann-Kendall Trend Analysis

In order to investigate the trend of cropland productivity and climate variables, the Mann-Kendall test (M-K test) [51], a non-parametric approach, was employed in this study. The M-K test has been widely used in vegetation trend analysis [52,53,54,55], and it is suitable for time series data of small sample size with occasional missing data and an assumption of normal distribution [51]. Two indictors—the M-K z-score and the median slope—were used. The z-score above or below 0 represents an increasing or decreasing trend, respectively. If there is a significant trend at the significant level of 0.05 (|z-score| ≥1.96), the rate of trend is then calculated by using the Theil-Sen’s slope [56]. More detailed information about the M-K test can be found in Neeti and Eastman [51].

2.6.2. Partial Correlation Analysis

Simple linear regression has been widely used for correlation analysis assuming dependent variables are independent of each other. However, climate factors are often found to be correlated with each other [2]. In order to detect how cropland productivity is driven by each climate factor at the temporal scale, a partial correlation analysis between climate variables and cropland productivity was performed. This analyzes the linear relationship between crop productivity and each climate factor after excluding the effect of other factors. A partial correlation coefficient (r) was calculated using Equation (3) to Equation (5) [57]. The value of r varies between −1 and 1 and shows the strength of the relationship between cropland productivity and a climate variable.

$$r_{xy}={{\displaystyle \sum }}_{i=1}^n\left[\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)\right]/\sqrt{{{\displaystyle \sum }}_{i=1}^n\left[{\left(x_i-\overline{x}\right)}^2{\left(y_i-\overline{y}\right)}^2\right]}$$

(3)

$$r_{xy,z}=\left(r_{xy}-r_{xz}{r}_{yz}\right)/\sqrt{\left(1-r_{xz}^2\right)\left(1-r_{yz}^2\right)}$$

(4)

$$r_{xy,zq}=\left(r_{xy,z}-r_{xq,z}{r}_{yq,z}\right)/\sqrt{\left(1-r_{xq,z}^2\right)\left(1-r_{yq,z}^2\right)}$$

(5)

$$t=r_{xy,za}/\sqrt{\left(1-r_{xy,za}^2\right)/\left(n-m-1\right)}$$

(6)

where, x, y, z, and a are the independent factors; r_xy is the zero order coefficient calculated using a simple linear regression between two variables (x and y); r_xy,z is the first order coefficient of the linear relationship between x and y after excluding the effect of z on x and y; r_xy,zq is the second order co-efficient of the linear relationship between x and y after excluding the effect of z and q on x and y, respectively.

The significance of the partial correlation coefficient was also evaluated with the Student’s t-test statistic (Equation (4)). A correlation was presented at three significance levels 0.01 (p-value, p = 0.01), 0.05 (p = 0.05) and 0.1 (p = 0.1), with decreasing strength of correlation, respectively.

3. Results

3.1. Assessment of CFAPAR for Cropland Productivity

To evaluate how well CFAPAR represents cropland productivity, it was necessary to investigate its relationship with cropland productivity. Figure 2 shows a reasonably good linear correlation between crop NPP (NPP_CAR) and CFAPAR at the CAR level (CFAPAR_CAR) on the Canadian Prairies. The correlation coefficient (r) is 0.60 (p < 0.01, n = 513). The result indicated that CFAPAR could be considered as a reasonable surrogate of cropland productivity, in agreement with some other studies [4,58].

3.2. Spatial Patterns of Cropland Productivity and Climate Indices

As discussed above, CFAPAR_CAR was used as a surrogate for cropland productivity to assess the variation of cropland productivity during 2000–2013. The spatial patterns of average annual CFAPAR from 2000 to 2013 on the Canadian Prairies are shown in Figure 3. In general, the spatial patterns of CFAPAR were in conformity with the distributions of agro-climatic zones (Figure 1b). The average annual CFAPAR in the sub-humid region is the highest among the three agro-climatic zones (~75), followed by the semi-arid region (~56) and the arid region (~50) (Figure 3). At the provincial level, MB had the highest CFAPAR (~85) among the three provinces, followed by AB (~80) and SK (~65).

The spatial distributions of climate indices, including cumulative rainfall, daily mean temperature and cumulative daily radiation during a growing season, are shown in Figure 4. In general, their spatial distribution patterns are consistent with the distribution of agro-climatic zones. The growing season cumulative rainfall in the prairie region was generally less than 450 mm (Figure 4a), and decreases from the sub-humid area to the arid area. The sub-humid area received the highest amount of rainfall in each growing season (~273 mm), especially in the area along the western and eastern borders of the region (>350 mm). The cumulative rainfall was higher in the semi-arid area (~205 mm) than in the arid zone (~183 mm). The growing season mean temperature generally decreased from the sub-humid to the arid region (Figure 4b). The pattern of temperature distribution in the sub-humid region was more complex than in the other two regions as the temperature in AB was lower than that in SK and MB. The cumulative radiation during a growing season showed an opposite order of distribution compared with that of daily average temperature (Figure 4c). The highest value of radiation occurred in the western portion of the sub-humid region in the province of AB.

By comparing the spatial patterns of cropland productivity in Figure 3 and climate indices in Figure 4, we observed that the spatial patterns of the climate variables over the growing season were very similar to that of cropland productivity. One of the factors that contributed to these spatial patterns was the duration of the growing season length, which was determined by SOS and EOS, as discussed by Meroni et al. [4]. The overall correlation between average annual CFAPAR and climate variables are shown in Figure 5. Both cumulative rainfall (r = 0.79, p < 0.01) and cumulative radiation (r = 0.76, p < 0.01) were positively correlated to cropland productivity, while the growing season mean temperature showed a negative relationship (r = −0.73, p < 0.01) with cropland productivity. Partial correlation analysis (

Table 1. Partial Correlation Matrix, a partial correlation analysis was conducted between average annual CFAPAR and average annual cumulative rainfall, average annual growing season mean temperature and average annual cumulative short wave radiation during 2000–2013. The dataset is shown in Figure 5.

	CFAPAR	Rainfall	Temperature	Radiation
CFAPAR	−1.00	0.60	−0.46	0.00
Rainfall		−1.00	0.39	0.49
Temperature			−1.00	−0.64
Radiation				−1.00

) demonstrated that spatial patterns of the annual CFAPAR was mostly affected by the cumulative rainfall (r = 0.60, p < 0.01) and the growing season mean temperature (r = −0.46, p < 0.05), and not by cumulative radiation (r = 0.00, p > 0.1).

3.3. Trend of Cropland Productivity

To examine the spatial-temporal variation of cropland productivity, we conducted the trend analysis using the M-K test at the pixel level. Figure 6a shows the trend of CFAPAR from 2000 through 2013 represented by the M-K z-score. Results show that the z-score in most crop areas is above zero (Figure 6a), indicating a general increasing trend of cropland productivity from 2000 to 2013. The z-score in some areas such as the southern and central regions of SK and AB was higher than in other cropland areas. Some small areas such as the area in the northern part of MB and northwestern part of AB, show a decreasing trend (z-score < 0). At the significance level of 0.05 (|z-score| ≥ 1.96), a significant increasing trend of cropland productivity was observed in AB and SK, as seen from the Theil-Sen median slope (>0.10) in Figure 6b. A significant decreasing trend of cropland productivity at the 0.05 level was observed in small areas scattered throughout MB (TS slope < −0.05).

To further analyze the trend of annual cropland productivity at larger scales, the M-K test was also performed at the provincial level and for agro-climatic zones. Figure 7 shows the trends of cropland productivity in the three provinces. Cropland productivity showed a strong increasing trend in AB (Z-score = 2.08 > 1.96, slope = 0.91 year⁻¹, Figure 7a) and SK (Z-score = 2.85 > 1.96, slope = 0.96 year⁻¹, Figure 7b). No significant trend was observed for cropland productivity in MB at the 0.05 significance level. Figure 8 shows that cropland productivity had a significant increasing trend in the sub-humid zone (Z-score = 2.74 > 1.96, slope = 0.71 year⁻¹), the semi-arid zone (Z-score = 2.63 > 1.96, slope = 1.23 year⁻¹) and the arid zone (Z-score = 2.52 > 1.96, slope = 1.20 year⁻¹), with the largest increase in the semi-arid zone and the smallest in the semi-humid zone. This indicates a general increasing trend of cropland productivity in the Canadian Prairies. Figure 8 also shows that cropland productivity in the sub-humid zone is always higher than that in the other two zones. By comparing Figure 7 and Figure 8, it is apparent that trends in cropland productivity were more clearly observed when evaluated using the climatological districts than using the administrative boundaries.

3.4. The Relationships between Cropland Productivity and Climatic Variables

The spatial patterns of average annual CFAPAR from 2000 to 2013 in Section 3.2 highlighted the potential connections between cropland productivity and climate variability. Figure 9 shows the partial correlation coefficients of the three climatic variables with annual cropland productivity from 2000 to 2013.

For cumulative rainfall in a growing season (Figure 9a), the variability of cropland productivity is generally sensitive to rainfall variability with a positive relationship in more than 70% of the study areas (Figure 9a). The partial correlation coefficient (r) over most cropland areas in the Canadian Prairies was higher than 0.55, especially in AB and SK. Cropland productivity in most areas was negatively affected by growing season mean temperature (Figure 9b); however, except for areas in the arid zone located mainly in southern SK and southeast AB, the effect was not significant. For the cumulative radiation during the growing season (Figure 9c), there was no significant correlation between cropland productivity and radiation at the 0.05 level in most areas, indicating that cropland productivity in these areas was not normally limited by radiation. However, cropland productivity in some localized croplands of central AB and northeast SK was sensitive to radiation with a positive partial correlation coefficient greater than 0.55 (p < 0.01, n = 14). For AB and SK, it could be observed that rainfall played the most important role in variation of cropland productivity with a positive relationship, followed by temperature variability and radiation variability (Figure 9). The variation of cropland productivity in MB was more affected by the temperature variability (p < 0.10) and the limiting factors were more complicated than that in AB and SK. In summary, the analysis in this section demonstrated that cropland productivity variability in most areas of the Canadian Prairies was mainly driven by the amount of rainfall.

To further determine the trend of cumulative rainfall during observed years, the analysis was conducted using the M-K test and the results are shown in Figure 10. Areas with a positive z-score, indicating an increasing trend, were distributed in most areas of the Canadian Prairie region. However, the increasing trend was not significant at the level of 0.05 (|Z-score| ≤1.96) except in some localized areas.

To further evaluate the relationship between cropland productivity and rainfall, the spatial distribution of the coefficient of variation (CV) of rainfall is shown in Figure 11 and compared with that of cropland productivity. Although there was no apparent trend for rainfall, the CV of rainfall showed strong fluctuation. Rainfall variability ranged from 40% to 60% across the three agro-climatic zones, with a higher variability in the semi-arid zone than in the sub-humid zone. The CV distribution of cropland productivity (Figure 12) shows a consistent pattern as in Chipanshi, et al. [12], and has a similar spatial pattern to that of rainfall (Figure 11). A strong and positive linear correlation between the CV of CFAPAR and that of rainfall (r = 0.66, n = 593, p < 0.01) was observed (Figure 13). Hence, these results indicate that the cropland productivity variability was partly determined by the rainfall variability. These results also further explain why the cropland productivity showed positive correlation to rainfall variability in most areas of the Canadian Prairies.

4. Discussion

Our results show that CFAPAR, a seasonal indictor of green area formation [38], has a strong linear correlation with the cropland NPP (r = 0.60, p < 0.01, n = 513). This is consistent with literature studies [4,32,38,59]. For example, Meroni et al. [4] found that the correlation coefficient (r) between the measured herbaceous dry biomass in Senegal and the CFAPAR retrieved from the SPOT-VEGETATION was 0.58, and López-Lozano, et al. [38] found that CFAPAR had a stronger correlation with official yields under water-limited conditions than that of non-water-limited areas. The main reason is that FAPAR plays an important role in photosynthetic process [30,31,38]. Thus, CFAPAR, which has a much higher spatial resolution than that of crop yield statistics, could be used to demonstrate the influence of climate variables on cropland productivity at different spatial-temporal scales. However, a larger portion of the NPP variability was accounted for by the accuracy of CFAPAR. The contribution factors to CFAPAR variability may include, but are not limited to, (a) the uncertainties in MODIS FAPAR product caused by mixed pixels and retrieval algorithm [32,42,60]; (b) a general cropland mask obtained from the MODIS land cover product [43]; (c) the uncertainties in the extracted phenological metrics [4]; and the modeling errors in converting census crop yield to NPP [37,50], especially the former two factors.

The results in this study also showed that cropland productivity had an increasing trend, especially in AB and SK (Figure 6 and Figure 7). This increasing trend has been attributed to the contribution of technology improvement [12]. Furthermore, the influence of cropland variability on cropland productivity was negligible in this study as the crop area data from Statistics Canada (http://www.statcan.gc.ca/start-debut-eng.html) showed that the variability of total cropland area on the Canadian Prairies was quite stable from 2000 to 2013. For the cropland productivity variation analysis, cropland productivity located in the semi-arid and arid zones showed strong variation (Figure 12). Our results are consistent with Chipanshi, et al. [12] who analyzed the variation of the surveyed yield from Statistics Canada over 1985–2012 at the CAR level. Further analyses about the relationships between CFAPAR and climate variables suggest that rainfall has played a key role in the interannual spatial-temporal variation of cropland productivity. Rainfall had a smaller effect on cropland productivity variation in MB than in AB and SK (Figure 9 and Figure 11). A major reason is that annual cumulative rainfall in MB is usually ample (Figure 4a) and its variation is smaller (Figure 12) than that in the other two provinces. The analysis of the spatial relationships between cropland productivity and rainfall suggest that rainfall could account for about 50% of the cropland productivity variation. Thus, it can be concluded that higher cropland productivity variation in the Canadian Prairies is often accompanied by higher variability in rainfall. Cropland productivity was also observed to be negatively affected by temperature in some areas of the Canadian Prairies, which is consistent with results of other researchers, who have indicated that water stress is one of the major yield limiting factors in the Canadian Prairies [2,3,6,8,15]. This also further confirms why the CFAPAR in the Canadian Prairies in the years of 2001, 2003, 2006, 2009 and 2012 were below the corresponding trend line (Figure 7 and Figure 8). Specifically, following the severe drought in 2001 (Fang et al., 2014; Hanesiak et al. 2011), cropland productivity in AB and SK (Figure 7a,b) dropped to the lowest level during the observed period. Cropland productivity in these areas was continuously influenced by the drought in the following two years and then began to recover in 2004. In 2009, Canada experienced the driest year of the past 70, and comparable to the drought year of 2002 (Canada’s Top Ten Weather Stories for 2009, https://ec.gc.ca/meteo-weather/), most cropland in central and eastern Canada was severely affected by the drought and cropland productivity in AB, SK, and MB was below the normal level of yield. In 2012, crop conditions were severely affected by a drought again and cropland productivity was below average (Canada’s Top Ten Weather Stories for 2012, http://www.ec.gc.ca/meteo-weather/).

Although the results from this study are encouraging, several uncertainty factors other than those investigated in this study are worth noting. The first uncertainty is from the MODIS land cover product, which has an accuracy of about 80% in identifying croplands globally [43]. SOS and EOS extracted from MODIS product may be affected by pixel mixing effect. This is inevitable when using coarse resolution MODIS products; however, the impact should be limited to the Canadian Prairies agricultural region as it is cropping intensive and major crops have roughly a similar phenology cycle. At the same time, the identification of crop SOS and EOS from remote sensing may be different from dates identified from field measurement or crop model simulation [61]. Therefore, the uncertainties in the assessment of the trend of cropland productivity and the influences of climate variables on cropland productivity in this area may also be affected by the length of growing season derived from MODIS FAPAR.

It should be acknowledged that our analyses were focused solely on the influence of cumulative climatic conditions, such as cumulative daily rainfall during the growing season, on cropland productivity variation. Therefore, our study provides an annual synoptic view of the climatic impact on cropland productivity in the Canadian Prairies. In reality, cropland productivity is impacted by climatic condition in a more complex manner and finer spatial-temporal scale than that used in this study. Climatic conditions may have different impacts on the final yield at different growth stages, and the impacts may also be significantly different in different regions [3,19,62]. For instance, spring wheat could benefit from a slight water stress at the heading-soft dough stage in Manitoba [2], whereas canola yield in Saskatchewan could be negatively impacted by low rainfall and high temperatures at the early stage of flowering [62]. Further studies assessing the effect of climate variables on crop growth stage and final yield at finer spatial scale using remote sensing are needed.

5. Conclusions

In this study, spatial patterns and temporal trends of cropland productivity and their relationships with climate variables from 2000 to 2013 in the Canadian Prairies were analyzed based on time series of MODIS-derived FAPAR and gridded climate datasets. The results effectively revealed the effects of climate variability on cropland productivity variation based on remote sensing observations in the Canadian Prairies with a better spatial coverage than using traditional in-situ data. The main conclusions can be summarized as follows:

(1)

The cumulative FAPAR during the growing season can be presented as a proxy of cropland productivity;

(2)

There was, in general, an increasing trend in cropland productivity during the period from 2000 to 2013 over most of the cropland area of the Canadian Prairies.

(3)

Temporal and spatial variabilities in cropland productivity are both connected to rainfall variability, with temperature being a negative factor in arid regions. The trend towards increasing cropland productivity was somewhat greater in the more arid regions of the Canadian Prairies.

This technique shows how FAPAR can be used to evaluate productivity trends on large spatial scales in response to climate and other stressors for crop production, in particular in areas where sparse information is available from other data sources on productivity or yield.