4.1.1 Precipitation and temperature data and the CDI

Gridded daily rainfall data from 1901–2004 available at 1^∘ × 1^∘ spatial resolution from the India Meteorological Department (Rajeevan et al., 2006) and gridded daily temperature data from 1969–2005 available at the same spatial resolution from India Meteorological Department are used in this study. Since the daily temperature data are available only for 37 years, we used the daily climatology, i.e., the mean daily temperature, for the remaining 77 years (Devineni et al., 2013). The daily climate time series grids were spatially averaged over the Satara district. This process resulted in a time series of daily precipitation and temperature estimates for 104 years. The daily reference crop evapotranspiration (ET₀) was developed based on the daily time series of minimum, mean and maximum temperature data and extraterrestrial solar radiation (Hargreaves and Samani, 1982). The Hargreaves method is used globally to predict ET₀ in regions where data availability is limited to air temperature data (Allen et al., 1998). Seasonal daily rainfall data from 2005 to 2013 for the Satara district were collected separately from a website maintained by the Agricultural Department of Maharashtra State and used to augment the 104 years of rainfall and temperature data. The CDI was computed for each of these 113 seasons using the daily rainfall data and reference crop evapotranspiration. This will serve as the predictand for our forecast model. We remind the reader that Fig. 1 illustrates the computation of CDI.

CDI as a water stress measure is a proxy of not only crop water stress but also irrigation and water storage requirements. Consider Fig. 1. When daily seasonal rainfall is low or when rainfall enters an inactive phase for a considerable period of time, as displayed by the vertical cyan bars, the amount of daily accumulated water deficit increases to reflect the disparity between water supplied as rainfall and the water required by the crop to sustain itself, as displayed by the red curve in Fig. 1. The highest point, or peak, on the black deficit time series in Fig. 1 is the value of CDI, and it prepares us for the worst-case scenario of a deficient water supply for the crop. This can be calculated for multiple crops, with each CDI value depending on the specific crop's water demand and the location and time of planting. This gives the stakeholder a conservative estimate of how much additional water is needed beyond what nature is willing to supply in order to maintain critical yields while apportioning water resources intelligently. Since agriculture tends to be one of the largest consumers of water – about 70 % of all the world's freshwater withdrawals go towards irrigation use (USGS, 2017) in addition to what is rainfed – this is an integral part of water resource management.

The annual time series of the CDI computed for the JJAS season (referred to as the Kharif season on the Indian subcontinent) in Satara is presented in Fig. 3. We have standardized the CDI values as the percentage difference each year from the 113-year average of CDI. The long-term average CDI for growing potatoes in Satara is 241 mm. This is equivalent to approximately 975.3 m³ of water used for irrigating a 4046.86 m² farm of potatoes on average throughout the season. The percent differences in Fig. 3 refer to percentages of this number, i.e., a 10 % increase in CDI indicates an additional requirement of 97.5 m³. From Fig. 3, it is clear that (a) Satara experiences recurrent droughts with intermediate wet periods and (b) there is year-to-year persistence in the incidence of these droughts. Such variations and epochal changes are typically modulated through large-scale global climate patterns. Investigating the relationship between the monsoon deficit and the large-scale climate teleconnections could enable the development of models that can be used to understand and predict the variability in the CDI in the region.

4.1.2 Climate precursors and climate data

Our goal was to develop a simple statistical model for predicting the CDI for potatoes grown in Satara. The generalized climate forecast models available at low spatial resolution are not specific enough for this task. Consequently, the first objective was to identify appropriate climate predictors before the monsoon starts in June. There is an extensive history of developing long-range predictions of monsoon rainfall that are based on various regional to large-scale climate predictors (Walker, 1924; Thapliyal, 1987). A variety of seasonal forecasts of the Indian summer monsoon rainfall (ISMR) are documented and available for reference (Gadgil et al., 2007; Kumar et al., 1995).

It is well established that inter-annual climate modes such as ENSO associated with anomalous sea surface temperature (SST) conditions in the tropical Pacific Ocean influence the inter-annual variability of ISMR (Parthasarathy and Pant, 1985; Shukla and Paolino, 1983). Anomalously warm tropical eastern Pacific SSTs (El Niño) are associated with a drier-than-normal ISMR, whereas anomalously cool tropical eastern Pacific SSTs (La Niña) are associated with a wetter-than-normal ISMR (Sikka, 1980; Parthasarathy and Pant, 1985; Rasmusson and Carpenter, 1983). Ihara et al. (2007) have suggested that the ENSO warm (cool) phases shift the location of the tropical Walker circulation and cause deficient (excessive) rainfall by suppressing (enhancing) the convection over India. Hence, ENSO indices were chosen to be among the candidate predictors for the forecast model. Raw monthly SST data for the Niño 3, Niño 4, Niño 12 and Niño 34 indices were taken from the Royal Netherlands Meteorological Institute (KNMI) climate explorer database (KNMI, 2014).

For each given raw ENSO index (3, 4, 12 and 34), we considered three different types of derived ENSO indices: a December–January–February (DJF) seasonal average, a March–April–May (MAM) seasonal average, and a MAM minus DJF (MAM − DJF) differenced time series. Among the Niño indices calculated, the change in the tropical Pacific SSTs from December to May (MAM − DJF trend) was found to be of significance by previous investigators. Shukla and Paolino (1983) found the correlation coefficient between the MAM–DJF trend pressure anomalies and the ISMR to be a significant −0.42. Their investigation showed that the Darwin pressure anomalies decrease from DJF to MAM before the occurrence of heavy monsoon rainfall and increase prior to the occurrence of deficit monsoon rainfall. Parthasarathy et al. (1988) found the correlation coefficient between this winter-to-spring trend and ISMR over the period 1951–1980 to be between 0.40 and 0.52 in magnitude, depending on the specific region within the tropical pacific. Hence, MAM-DJF trends from Niño 3, Niño 4, Niño 12 and Niño 34 were considered to be potential model predictors. Parthasarathy et al. (1988) found that the MAM-averaged tropical Pacific SSTs over the box 14 to 20^∘ N, 176^∘ E to 160^∘ W had a correlation of −0.40 with ISMR, convincing us to consider this average as well. In addition to the MAM and MAM − DJF averages, we computed the winter season (DJF) average, although DJF-averaged tropical Pacific SSTs were not found to be significant in the literature. However, it is worth noting that Parthasarathy et al. (1988) found that the correlation coefficient between the Darwin SLP during the DJF season and ISMR was +0.39. As the concurrent season (JJAS) state of ENSO has an important, well-documented impact on ISMR, we also elected to include the Niño 34 JJAS average. As mentioned earlier, an El Niño event during the JJAS season is strongly associated with an anomalously dry JJAS rainfall season in India, while a La Niña event during the JJAS season is strongly associated with an anomalously wet JJAS rainfall season in India, prompting our choice. We coupled the JJAS seasonal average for the Niño 34 index with forecasts of the JJA and JAS seasonal averages for the Niño 34 index. These forecasts were obtained from the International Research Institute for Climate and Society (IRI) ENSO forecast page and covered the period 2002–2013. These forecasts can be used to forecast JJAS monsoon CDI in place of the observed Niño 34 JJAS values on a real-time basis. These forecasted values were averages of the projections from at least six distinct statistical/dynamical models, with one average for the JJA season and one average for the JAS season. Together, we start with a total of 13 ENSO-based indices.

Figure 4Spearman's rank correlation between the CDI in Satara and SST field during the same JJAS season. SST region in the Indian Ocean (red box) that influences the CDI has a statistically significant correlation at the 95 % significance level.

Download

Other candidate predictor variables include concurrent season (JJAS) eastern Indian Ocean SSTs known as the Indonesian Throughflow or ITF. Warm, low-salinity water from the Pacific is introduced into the Indian Ocean via the ITF and is considered to be an integral component in the heat and hydrological budget of the Indian Ocean (Gordon et al., 1997). The ITF waters are also believed to influence SSTs and associated ocean–atmosphere coupling within the Indian Ocean, making it an important aspect of monsoon climate research (Gordon et al., 1997). Thus, the ITF was also selected to be a candidate predictor in the model. During the JJAS monsoon season, the ITF is strengthened considerably, allowing an abundant amount of relatively warm water to be injected into the Indian Ocean. Eastern Indian Ocean SSTs during the JJAS season correspond to enhanced (suppressed) atmospheric convection during the anomalous warming (cooling) of the Indian Ocean waters, which in turn supplies (robs) the developing monsoon of much-needed moisture. We found that the Spearman's rank correlation coefficient between CDI in Satara and the average SST anomalies over 20^∘ N and 5^∘ S and 100 and 130^∘ E (the region representing ITF) during the JJAS season is around −0.35 (statistically significant at the 95 % level), suggesting that warm conditions in the ITF region result in below-normal CDI, or low crop water stress. Figure 4 presents the field correlation map of SST anomalies with CDI. For these reasons, we chose the concurrent season ITF data to be a candidate predictor. The ITF data were collected from the IRI data library and consist of two components, namely an observed component and a forecasted component. The observations consist of measured eastern Indian Ocean SST anomalies during the JJAS season from 1901 to 2013. The forecasts consist of JJAS-season ITF values retrospective of the ECHAM4.5 global climate model and cover the period 2001–2013. Skillful forecasts for the tropical SSTs based on coupled ocean–atmosphere general circulation models have been in operation from various climate centers since 1998. Hence, in the forecasting scheme, we used the ITF derived from the forecasted SST state issued in May from the ECHAM4.5 operational forecasting center (available from IRI data library: http://iridl.ldeo.columbia.edu/SOURCES/.IRI/.FD/.ECHAM4p5/.Forecast/.ca_sst/.ensemble24/ (last access: 5 February 2017); Li and Goddard, 2005; van den Dool, 2007; Roeckner et al., 1996). The observed JJAS ITF data are used to train the model, while the retrospective JJAS ITF forecasts are used to make forecasts for the years 2001–2013.

4.2.1 Predictor selection

Given a pool of candidate predictors, the next step is to select the best subset of those predictors. The predictors used in the forecasting model were chosen based on an exhaustive search method. In the exhaustive search method, all possible combinations of the candidate predictor variables are used to develop models that are cross-validated on historical data. Skill metrics are then used to compare the predictive accuracy of each combination. In the present study, we began with 113 years of CDI data and fourteen candidates: Niño 3 DJF, Niño 3 MAM, Niño 3 MAM − DJF, Niño 4 DJF, Niño 4 MAM, Niño 4 MAM − DJF, Niño 12 DJF, Niño 12 MAM, Niño 12 MAM − DJF, Niño 34 DJF, Niño 34 MAM, Niño 34 MAM − DJF, Niño 34 JJAS and ITF. The exhaustive search method utilized the k-NN cross-validation algorithm and 40 years of training data (1901–1940) to build forecast distributions for each of the years 1941–2013. At each step, the training data were updated to include data from all of the years up to the year being cross-validated. Thus, we always only use the historical data and update the model each year with the information from the previous year, much as a regular user of the forecast system would have to do. These forecasting distributions, built over a 73-year record (1941 to 2013) were created successively for every unique combination of two variables, every unique combination of three variables, and so on until we reached the entire pool of predictors.

For each and every possible unique combination of the predictor variables, we obtain a matrix of 73 columns. For each of these 73 years, the squared error and ranked probability score (Epstein, 1969; Murphy, 1969, 1971; Candille and Talagrand, 2005) were computed, and from this the root-mean-squared error (RMSE) and ranked probability skill score (RPSS) were computed. In this manner, a single RPSS value and RMSE value were calculated for every possible combination of the predictor variables. We chose the following combination of predictors based on the relative optimality of both their RPSS and RMSE scores: Niño 12 MAM − DJF, Niño 34 MAM − DJF and ITF, and this set of variables had an RMSE of 49.25 mm of required (JJAS) seasonal water storage and an RPSS of 0.26. We devised a simple but effective decision rule for determining the optimal choice of predictors based on ranking the metric values. This is especially useful when the number of combinations of variables is unwieldy. Optimality was determined by assigning a rank number to the RMSE and RPSS values in such a way that the first number was assigned to the lowest RMSE value, the second to the second lowest RMSE value, and so on, and the first number was assigned to the largest RPSS value, the second to the second largest RPSS value, and so on. For a fixed number of cross-validated predictor candidates for each RMSE/RPSS pair and for one pair of each combination of predictors, we determined an RMSE and RPSS rank and took the sum of these ranks. The smallest of these sums corresponds to the best or optimal set of predictors among all possible sets of cross-validated predictors. We then compared the ranked sum while considering the number of predictors in order to choose the best set of predictors. The chosen trio of predictors mentioned above had the unequivocally highest value of RPSS and second lowest RMSE value out of all possible combinations of the original set of 17 candidates, the lowest RMSE being only slightly smaller at 48.92 mm. Conceptually, this procedure is similar to the “best subsets regression” or “step-wise regression” (Helsel and Hirsch, 2002), but in the spirit of using the k-NN algorithm for forecasting, we designed this selection scheme to use the k-NN algorithm instead.

CDI forecasts were subsequently made using the selected set of predictors. The forecast procedure is tested using the leave-one-out cross-validation method. Each historical observation is omitted in turn, and the model is developed using the remaining years of data. A prediction of the observation that was not kept in the model-building set is then made and compared with the actual outcome for that year. Results from a variant of this approach are presented in the next section. The CDI for the 2001 Kharif season is predicted using the model developed based on data from 1901 to 2000. Similarly, the CDI for 2002 is predicted based on the model that is developed using the data from 1901–2001. Thus, as we move from year to year, we update the model observations and predict the future state.

4.2.2 The k-nearest neighbors real-time forecasting model

The forecasts were developed using a non-parametric k-nearest neighbors (k-NN) model. This is a data-driven approach that develops a conditional probability distribution of the CDI given the predictors by first identifying the k-historical climate conditions that are most similar to the current values of the climate predictors and then randomly drawing the vector of CDI values in the historical data that correspond to these k neighbors. The neighbors are weighted so that the closer or more similar neighbors are chosen more often than those further away. The key steps are as follows.

Let X be the design matrix of size n×p, where p = number of predictors selected from the original pool of candidates. Let x_i denote the ith row of X. Hence, x_i is a vector containing the values of each of the p predictor variables during year i. In denoting the current values of the predictors by x_c, the idea is to find k such predictor vectors from the historical record (i.e., find k values of x_i with i<c) that are most “similar” to the value of x_c and use this information to construct a sampling distribution of CDI from which we can issue probabilistic forecasts. The number of neighbors in the model, or k, represents the number of degrees of freedom in the model, and should be chosen with care, as the choice of k affects the skewness and level of uncertainty in the sampling distributions. After trying several different values for k, we found an optimal value to be k=25. Rajagopalan and Lall (1999) recommend that, as a rule of thumb based on asymptotic arguments, k be roughly equal to $\sqrt{n}$, where n = the total number of observations. In our situation, it was evident that we required more neighbors than this rule would allow due to the skewness and variance apparent in the sampling distributions when using only 11 or fewer neighbors. Lall and Sharma (1996) note that if their discrete kernel is used for resampling the conditional bootstrap, then the weights for further neighbors will decrease. Hence, choosing a larger k may reduce the variance in the estimate while potentially increasing the bias in the estimate of the conditional distribution. Cross-validation can also be used to choose an optimal value for k in a given setting.

Let y be the n-dimensional vector of seasonal CDI values, each component of which represents the aggregate water deficit level over the JJAS growing season of every year in the historical record. Assume that y has been centered and normalized by its historical average to produce mean-normalized anomalies. The first step was to consider the individual distance values (with a specified metric) between x_c and x_i for i = 1, …, c−1. The chosen distance metric for our k-NN model was the Mahalanobis distance (Mahalanobis, 1936), represented as

where Σ is the covariance matrix of the training values in X. The Mahalanobis distance measure judges point separations in a metric space based on statistical dissimilarity, as opposed to a solely physical distance. Hence, the level of similarity between predictor values across different years is determined by the orientation and location of each point relative to the scatterplot of the predictor data. Large distances from x_c represent predictor values that are statistically anomalous in the context of the predictor data. After the Mahalanobis distances had been calculated, the k-smallest distance values, with k=25, were selected and the corresponding years in which these distances occurred were noted. These years, hereby referred to as the analog years, are the years during which the predictor signals were the most similar to those of the current year. The vector-valued predictors during these analog years are referred to as the neighbors of x_c.

The final step was to resample CDI values from the analog years. The resampling technique employed is a nonparametric method known as the bootstrap (Efron, 1979; Efron and Tibshirani, 1993). The idea behind the bootstrap component is to sample with replacement from a pool of data using the underlying distribution that generated the data to guide the sampling process. We chose not to assign a parametric family of distributions to the CDI data and instead estimated its underlying distribution non-parametrically using a kernel density estimator. This non-parametric method of k-NN bootstrapping was first introduced in Lall and Sharma (1996). Applications of the methods using different variants have since been presented (for example, see Rajagopalan and Lall, 1999; Souza and Lall, 2003 and references therein). We employed the same discrete resampling kernel proposed in Lall and Sharma (1996), which has the general form $K\left(j\right)=\mathrm{1}/(j\cdot S)$ with $S=\sum _{j=\mathrm{1}}^k\mathrm{1}/j$, where j is the rank of each neighbor of x_c, a rank of j=1 is assigned to the closest neighbor and a rank of j=k is assigned to the most distant neighbor. Our strategy was to build this kernel density estimator based on the ranks of the selected neighbors and resample the predictands from these analog years. We resampled from the 25 analog CDI values 1000 times, and each of the 25 values was resampled proportionally to the probability of its occurrence as determined by the density estimator.

4.2.3 Analyzing the k-NN results

The way in which model results are interpreted and presented is important for potential stakeholders. In this case study, our interest was in forecasting the CDI for a given potato growing season in Satara. The information from these forecasts can be of great use to potato farmers in Satara as well as corporations with investments in these farming areas. This necessitates a clear and concise communication of the forecast results.

The output of the k-NN model was a time series for each forecasted year consisting of 1000 realizations. This is the sampling distribution for the CDI and consists of mean-normalized anomaly values from the analog years converted to percentage values. As stated in the previous section, the deficit value from each analog year in the sampling distribution is represented proportionally to its probability of occurrence as assigned by a kernel density estimator. The sampling distribution is used to issue one-season-ahead probabilistic forecasts (i.e., the likelihood of a deficit for the forthcoming growing season). There are a whole slew of possibilities when it comes to using these sampling distributions for probability-based forecasts. Our approach consists of the following steps for a given forecasted growing season:

1. A box plot depicting the sampling distribution with the observed percent anomaly value superimposed on the box plot for every growing season forecasted. In using predictand anomalies, the historical mean becomes the zero line in the coordinate plane of the box plot.

2. A three-category forecasting system with the categories “above normal”, “normal” and “below normal” is used, provided that the historical mean and climatology are the threshold that is desired.

3. The probabilities for the categories specified in step 2 from the sampling distribution generated in step 1 are calculated and used to evaluate the accuracy and strength of the forecast based on contingency metrics such as hit rates and false alarms.

4. To get a sense of the spread and variability in the box plot distribution, the Interquartile Range (IQR) is calculated.

5. The value of the observed percent anomaly of the predictand is compared with the category in which the majority of the probability mass of the sampling distribution lies. This is of central importance in getting a basic sense of the accuracy of the forecast.

In general, the construction of such a sampling distribution allows the investigator the freedom to calculate probabilities on many different thresholds. The thresholds should be defined by the particular application and the needs of any stakeholders involved.