Yield estimation and sowing date optimization based ... - Springer Link

Climatic Change (2010) 98:565–580 DOI 10.1007/s10584-009-9746-4

Yield estimation and sowing date optimization based on seasonal climate information in the three CLARIS sites T. d’Orgeval · Jean-Philippe Boulanger · M. J. Capalbo · E. Guevara · O. Penalba · S. Meira

Received: 23 June 2008 / Accepted: 3 August 2009 / Published online: 8 October 2009 © Springer Science + Business Media B.V. 2009

Abstract The present article is a contribution to the CLARIS WorkPackage “Climate and Agriculture”, and aims at testing whether it is possible to predict yields and optimal sowing dates using seasonal climate information at three sites (Pergamino, Marcos Juarez and Anguil) which are representative of different climate and soil conditions in Argentina. Considering that we focus on the use of climate information only, and that official long time yield series are not always reliable and often influenced by both climate and technology changes, we decided to build a dataset with yields simulated by the DSSAT (Decision Support System for Agrotechnology Transfer) crop model, already calibrated in the selected three sites and for the two crops of interest (maize and soybean). We simulated yields for three different sowing dates for each crop in each of the three sites. Also considering that seasonal forecasts have a higher skill when using the 3-month average precipitation and temperature forecasts, and that regional climate change scenarios present less uncertainty at similar temporal scales, we decided to focus our analysis on the use of quarterly precipitation and temperature averages, measured at the three sites during the crop cycle. This type of information is used as input (predictand) for nonlinear statistical methods (Multivariate Adaptive Regression Splines, MARS; and

T. d’Orgeval · O. Penalba Centro de Investigacion del Mar y la Atmosfera, Buenos Aires, Argentina T. d’Orgeval Laboratoire de Météorologie Dynamique, Paris, France J.-P. Boulanger (B) Laboratoire d’Océanographie et du Climat, Expérimentation et Analyse Numérique, Paris, France e-mail: [email protected] M. J. Capalbo · E. Guevara · S. Meira Instituto Nacional de Tecnología Agropecuaria, Pergamino, Argentina

566

Climatic Change (2010) 98:565–580

classification trees) in order to predict yields and their dependency to the chosen sowing date. MARS models show that the most valuable information to predict yield amplitude is the 3-month average precipitation around flowering. Classification trees are used to estimate whether climate information can be used to infer an optimal sowing date in order to optimize yields. In order to simplify the problem, we set a default sowing date (the most representative for the crop and the site) and compare the yield amplitudes between such a default date and possible alternative dates sometimes used by farmers. Above normal average temperatures at the beginning and the end of the crop cycle lead to respectively later and earlier optimal sowing. Using this classification, yields can be potentially improved by changing sowing date for maize but it is more limited for soybean. More generally, the sites and crops which have more variable yields are also the ones for which the proposed methodology is the most efficient. However, a full evaluation of the accuracy of seasonal forecasts should be the next step before confirming the reliability of this methodology under real conditions.

1 Introduction The present article is a contribution to the CLARIS WorkPackage “Climate and Agriculture”, and aims at testing whether it is possible to predict yields and optimal sowing dates using seasonal climate information for three sites (Pergamino, Marcos Juarez and Anguil), representative of different climate and soil conditions in Argentina. Recent developments of climate models led to significant improvements in seasonal 6-month precipitation and temperature forecasts. It was proven that ensemble simulations run by climate models, with current oceanic state as initial condition, have some skill in forecasting probabilistic precipitation and temperature anomalies for several months in different regions (Barnston et al. 2003). The potential use of seasonal forecasts in crop production has been an active area of research in recent years. Applications are not straight-forward, and the information provided may not be perceived as sufficiently accurate or relevant (Vogel and O’Brien 2006; Hartmann et al. 2002). An issue raised in different operational frameworks is the communication between end-users, planners—local or federal governments, institutions—and climate scientists (Lemos et al. 2002; Roncoli 2006). The way seasonal prediction is presented by scientists, disseminated by scientists and planners, and used by planners or farmers strongly affects the global impact that seasonal prediction has on agricultural outputs, leading sometimes to insignificant changes in the production, and to a negative perception by the end-users (Lemos et al. 2002). During the 1990s, seasonal forecasts were mainly limited to forecasts of the El Niño Southern Oscillation (ENSO) and its impacts on large scale climatic patterns. Given the recent improvements in climate models, seasonal forecasts are now relatively skilled in predicting quarterly temperature and precipitation averages, which are often reduced to terciles. The lack of skills of the seasonal forecast systems at higher frequency does not yet allow the direct use of mechanistic crop models, which require daily data of precipitation, minimum and maximum temperature and solar radiation. In order to avoid any negative perception and favor the adoption of “state-ofthe-art” seasonal forecasts, it is essential to estimate accurately their benefit for


567

crop production. In this article, we focus on three pilot sites located in Argentina: Pergamino (Province of Buenos Aires), Marcos Juarez (Provinde of Cordoba) and Anguil (Province of La Pampa). These three pilot sites were selected because of contrasting present climates and different responses to the 1970s observed recent climate change. Atmospheric circulation over these regions is known to be influenced by ENSO (Vargas et al. 1999; Grimm et al. 2000; Ropelewski and Halpert 1996). Recent studies of the impact of climate variability on seasonal climate prediction in this region therefore mainly focused on ENSO impacts (ENSO and crop yields, Podestá et al. 1999; ENSO and land allocation, Messina et al. 1999). However, using ENSO seasonal prediction only for agricultural decision-making might prove limited in this region. Letson et al. (2001) conducted farmer surveys in Pergamino. They found that, in addition to the doubts about the reliability of the forecast, the farmers’ “incomplete” knowledge of how ENSO affects their region may also pose an obstacle to further use of climate information”. Bert et al. (2006) carefully detailed the decision-making process related to climate information in maize crop production in the Argentinean Pampas. Though some limits inherent to the study were indicated, they found that advice given by agricultural experts based on the predicted ENSO phase did not optimize the farmers’ margin. In fact, their retrospective simulations showed that no significant benefit could be gained by changing management options according to the ENSO phase. One main reason for the inefficiency of information based on ENSO is the wide range of climate variability during a given phase. This point is specifically addressed in a study of observed precipitation in the Paraná-La Plata hydrological basin by Boulanger et al. (2005), which showed that “the displacement of the teleconnection patterns from one [ENSO] event to another [impedes] the definition of a robust statistical relationship between ENSO and precipitation in the Paraná-La Plata basin”. They emphasized the need for more complex climate tools than linear statistical forecast systems for impact studies on society. These different results led us to the conclusion that an improved use of seasonal climate information for crop prediction in the region could be to rely on quarterly precipitation and temperature averages instead of ENSO phases. This article then focuses specifically on crop yields in three CLARIS pilot sites (Pergamino, Marcos Juarez and Anguil), representative of different climate and soil conditions in Argentina, different sensitivity to ENSO and different responses to the regional climate changes observed during the last 30 years (increase in precipitation leading to a westward displacement of cropping systems where climate conditions became more favorable). Marcos Juarez and Anguil, in particular, greatly benefited from increased precipitation, improving crop production. The aim of the article is to propose and test a methodology to turn seasonal climate information into a practical tool for agricultural decision-making before planting. The article does not address the accuracy of seasonal forecasts specifically. Its aim is rather to evaluate if typical well-predicted outputs can be turned into usable information for farmers in these agricultural regions. To model the relationship between climate and yields, nonlinear statistical methods were chosen because no grounds exist to suppose that the relationship is linear. Moreover parametric statistical methods were avoided to prevent the authors from subjectively choosing the shape of this relationship. To restrict the scope of the article, the present study only aims at evaluating the accuracy with which yields can be forecast and optimal sowing date can be chosen,

568


using quarterly temperature and precipitation averages combined with non-linear statistical methods in the three CLARIS sites. First, the estimated yields for each site and each crop provide valuable information that may be used by farmers to allocate land, select the type of seed (long vs. short cycles, resistance to specific conditions or threats,...), etc... Second, the optimal sowing date gives an indication for farmers on how much they can expect to improve yields using seasonal climate information. Section 2 describes the data, the model and the statistical methodology applied in this study. The use of seasonal climate information in yield forecasting is analyzed in Section 3. Section 4 analyzes the potential optimization of the sowing date using quarterly climate information. A discussion follows in Section 5.

2 Data, model and methodology 2.1 Data The three chosen sites for the CLARIS project correspond to locations for which the crop model has been calibrated and validated, and where agricultural management is best understood. Details on the sites are given in Table 1. The atmospheric data come from a dataset described by Penalba and Vargas (2004) with the corresponding checks and corrections. Radiation is estimated from sunshine duration observations Precipitation and temperature data are observed directly. There is almost no missing data and the dataset covers the 1959–2005 period for every site. Forty-six crop seasons can therefore be simulated. The three sites constitute a first overview of the different conditions faced by crop producers in an agricultural region covering the humid region of Pergamino and the drier and colder region of Anguil. Moreover, there are some differences in the soil types found on the respective sites. The soil characteristics are given by a soil survey map from INTA (see Table 1). The finer soil is found at Marcos Juarez with 6% of sand and the coarser at Anguil with 66.9% of sand at the surface. Typical crops for Pergamino and Marcos Juarez are soybean, maize and wheat. In Anguil, sunflower is also produced. Crop rotations are not detailed as they are not taken into account in this study. For each of the three sites, we decided to focus the study on the main summer crops (maize and soybean). 2.2 Crop simulations First, considering that our focus is on the use of climate information only, and that official long time series of yields are not always reliable (yields are often underestimated in order to pay less taxes) and are influenced by both climate and technology

Table 1 Characteristics of the three sites: Pergamino (PE), Marcos Juarez (MJ) and Anguil (AN) Site Coordinates Lat PE MJ AN

Lon

Mean climate

Surface characteristics (%)

P (mm d−1 –mm y−1 ) T (K) R (W m−2 ) Clay

33.56◦ S 60.33◦ W 2.6–950 32.42◦ S 62.09◦ W 2.4–875 36.30◦ S 63.59◦ W 2.0–730

289 291 288

16.3 16.0 15.2

22.7 25.1 14.2

Silt

Sand

64.8 68.9 18.9

12.5 6 66.9


569

changes, we decided to build a yield dataset with yields simulated by the DSSAT crop model (Decision Support System for Agrotechnology Transfer; Jones et al. 2003) already calibrated at the selected three sites and for the two crops under study (maize and soybean). Two specific crop models are used within DSSAT in order to simulate the two crops of interest: CERES for maize and CROPGRO for soybean. The model is used to simulate crop growth for a single crop at once without taking into account crop rotation, as the aim of the study is to focus on climate impact on two specific crops at three sites rather than on the optimization of crop rotation on a given site. For each crop and each of the three CLARIS sites, 46 years are simulated, using the daily forcing described in the previous section. The main parameters used for each genotype are the same as in previous validations of DSSAT for this region (Guevara and Meira 1995, 1999; Guevara et al. 1998, 1999; Meira et al. 1999). The soil state is initiated each year with soil wetness at 80% and nitrogen at 100% for every crop and every site. Then other management parameters are chosen according to in situ knowledge gained in the validation studies performed in each region. The main parameters are given in Table 2. The crops are not irrigated during the simulations because this is the usual practice in the Pampas. Moreover, simulations with irrigation have already been performed, and it is clear that climate impacts on crop yields are much lower when irrigation is available. The only varying parameter in the different simulations is the sowing date. For each site and each crop, an optimal sowing window is defined according to local knowledge. The optimal sowing window is 40 days long for maize and 30 days long for soybean at each site. In order to get different sowing options and to keep the analysis simple, three sowing dates are tested for each crop and each site. The first sowing date is the first day of the sowing window, the second is in the middle and the third is the last day. The tested sowing dates are therefore every 20 days for maize and every 15 days for soybean (see Table 2). 2.3 Statistical methodology 2.3.1 Climate information input The aim of the paper is to analyze the potential interest of seasonal climate information. The way this information is presented and used in the statistical analysis follows

Table 2 Management options used in the simulations with DSSAT at Pergamino (PE), Marcos Juarez (MJ) and Anguil (AN) for maize and soybean Crop

Site

Maize

PE MJ AN PE MJ AN

Soybean

a Default

Sowing dates 1

2

3

09/10b 09/10b 09/30 10/20 10/25a 10/30b

09/30 09/30a 10/20a 11/05b 11/10 11/15

10/20a 10/20 11/10b 11/20a 11/25b 11/30a

sowing dates (Sdef ) sowing dates (Salt )

b Alternative

Fertilization (kg N ha−1 )

Genotype

75

DK688 DK688 DK752 MG V MG IV MG III

0

Density (m−2 ) 7.5

30

570


the way it is presented in seasonal forecasts produced before the potential sowing dates, which is two-fold: –

Quarterly average of precipitation (P) and temperature (T) over the 6 months following September: • • • •

–

October–November–December (OND) November–December–January (NDJ) December–January–February (DJF) January–February–March (JFM)

Classification of these averages into terciles: below normal (−), normal (=), above normal (+).

A quarter is the typical time step for which seasonal forecasts show prediction skills (Barnston et al. 2003). Moreover, if seasonal forecasts are usually presented through classification into terciles, it makes sense to use directly the quarterly averages in our methodology for two reasons: (1) quarterly averages can easily be retrieved from national weather services and specialized research centers and (2) if one needs to estimate yields or at least if above-normal or below-normal yields are expected, it is best to build a statistical method to transform seasonal climate averages into yields and then to classify yields into terciles rather than do the opposite, as the transformation is expected to be non-linear. The statistical analyses described in Sections 2.3.3 and 2.3.4 and carried out in Sections 3 and 4 are also conducted using monthly values of Nino 3.4 index (Trenberth 1997) as climate input, without significant results. Therefore no relevant relationship between simulated crop yields or optimal sowing dates and ENSO information is found with our methodology. This is discussed in Section 5. 2.3.2 Default yield and optimal sowing date For each crop and each site, three simulations with different sowing dates are made. Based on these simulations, a default Sdef sowing date is defined as a constant sowing date for which the yields for the 46-year period are optimal. The default yield is the series of yields obtained with the sowing date Sdef . An alternative Salt sowing date can also be determined, as the best alternative to use on any given year when Sdef does not lead to optimal crop yields (indeed the optimal date for the 46-year period is certainly not the optimal date each year). Given the above, the aim of the article is to test whether quarterly precipitation and temperature averages combined with non-linear statistical methods can lead to: – –

Accurate forecasts of annual crop yields using Sdef as sowing date (default yields), and Optimizing the choice between Sdef and Salt for the sowing date in any given year.

This framework (optimizing the choice between a default sowing date and an alternative one) was selected because the methodology to choose between three sowing dates would be a lot more complicated than to choose between the two best ones, and would have, therefore, less operational interest. Indeed, the value of changing the sowing date does not only depend on the best sowing date, but also on the yield improvement that can be expected due to the change. A methodology is proposed here that takes quantitatively into account the yield improvement when


571

changing the sowing date, even if it does not allow choosing between more than two sowing dates (the default and the alternative—see Section 2.3.4 for an explanation of how yield improvement is taken into account). In order to select the “best” default and alternative sowing dates, a simple criterion based on the empirical percentiles of simulated yield is used. The “best” sowing date is defined to be the one that maximizes the yield corresponding to a certain percentile in yield. We optimized different percentiles such as 20% and 50%. Using the 20% percentile is useful for strongly risk-averse farmers (choosing preferentially to reduce losses during the worst years). Using the 50% percentile optimizes the mean/median over a long period of time. Both (whichever is useful to different farmer profiles) led to the same value of the default sowing date (see the value in Table 2 for each crop and each location). 2.3.3 Default crop yield forecast In order to predict crop yields (target variable) with quarterly precipitation and temperature averages (predictor variables or predictors), Multi-Adaptive Regression Splines—MARS (Friedman 1991)—are used. MARS is a nonparametric and nonlinear regression procedure using basis functions of the form max(x − t,0), max(t − x,0), min(x − t,0) or min(t − x,0) where x is a selected predictor variable, t is called a knot and represents a threshold above or below which the slope of the relationship between the predictor variable x and the target variable changes. Here MARS is used with no interaction between predictors (degree of interaction fixed to 1) which means that the models built are simple linear sums of basis functions, and are equivalent to piecewise multi-linear regression. In order to build the statistical models, 1960–1990 is used as training data. 1991– 2005 is used as validation data to test the robustness of the model. Various partitions between training and validation have been tried and this choice is found to be the best compromise to get a training dataset long enough to build MARS models and a sufficient validation dataset to test the robustness of models. Finally, a penalty of 3 is applied in model construction. The penalty indicates how much the algorithm penalizes the introduction of new basis functions in the model. A penalty of 3 with a degree of interaction 1 is a choice that allows a reasonable limitation of the number of predictor variables chosen by the algorithm. This is done to avoid over-fitting, which results in a good fit between the model and the training dataset, but no skill to predict the target variable outside this training dataset. 2.3.4 Optimal sowing date forecast Optimizing the sowing date also requires a non-linear method. However, the variable to predict is a binary one corresponding to the choice between Sdef and Salt . Therefore, regression-like methods such as MARS do not seem appropriate because the variable to predict is not continuous. Instead, recursive partitioning is used to build a classification tree. The goal of a classification tree is to predict the response of a categorical variable (here the choice between default and alternative sowing date) based on predictor variables. It is often considered that MARS is a form of linearization of classification trees (Friedman 1991) because a basis function in MARS uses a predictor to classify the sensitivity of the target variable (sensitive/insensitive on

572


each side of a knot). Examples of classification trees are provided in Breiman et al. (1984). Here, classifications of precipitation and temperature into terciles are used as predictor variables rather than precipitation and temperature themselves. This choice is consistent with the way seasonal forecasts are currently produced. The resulting analysis checks the accuracy of sowing decisions based on information such as “Above-normal (upper tercile) precipitation in October–November–December”. More examples will be given in Section 4.2 with the obtained classification trees. Similarly to what is done for MARS models, classification trees are built with 1960–1990 as training data and 1991–2005 as validation. Finally, in the training period, each evaluation of the variable to be predicted (choice between two sowing dates) is weighted by the difference in yield between both sowing dates. The aim is to put more weight on picking the right sowing date when the resulting yield improvement is greater.

3 Yield forecast 3.1 Analysis of yield variability First, simulated default yields (corresponding to the Sdef sowing date) for the different crops and locations are plotted in Fig. 1. The yields at Anguil are clearly lower than at Pergamino and Marcos Juarez, especially for maize: the average yield

1970 1980 1990 2000

Yield (t/ha) 1970 1980 1990 2000

1960

Year

25 15

Yield (t/ha)

15

5 0

0

5

10

Yield (t/ha)

20

25

f. Soybean - Anguil

20

25

e. Soybean - MarcosJuarez

20 15 10 5 0

1970 1980 1990 2000

1970 1980 1990 2000

Year

Year

d. Soybean - Pergamino

1960

10 20 30 40 50 60 70 0

1960

Year

10

1960

Yield (t/ha)

c. Maize - Anguil

0

Yield (t/ha)

10 20 30 40 50 60 70

b. Maize - MarcosJuarez

0

Yield (t/ha)

10 20 30 40 50 60 70

a. Maize - Pergamino

1960

1970 1980 1990 2000

Year

1960

1970 1980 1990 2000

Year

Fig. 1 Time series of reference yields (in tons per hectare) from the DSSAT model for maize and soybean at Pergamino, Marcos Juarez and Anguil between 1960 and 2005


573

at Anguil is nearly two thirds of the one simulated at the other two sites. Yields at Anguil are also the most variable ones: the standard yield deviation is 40% above of the average yield for both crops. This implies that yields can easily change by a factor 5 from 1 year to another. At Pergamino and Marcos Juarez, maize yields are less variable than soybean yields. For both crops, the average yields are of the same order for Pergamino and Marcos Juarez. Maize yields are clearly less variable at Marcos Juarez whereas soybean yields are only slightly less variable at Pergamino. However, very low yields are simulated even for maize at Marcos Juarez (Y = 0 in 1994) and soybean at Pergamino. Therefore yield variability is fairly high for every site and every crop. The next subsection determines which part of this variability can be predicted with quarterly temperature and precipitation averages. 3.2 Statistical forecast of yields

1960 1970 1980 1990 2000

40 50 60 70

Yield (t/ha)

0 1960 1970 1980 1990 2000

1960 1970 1980

1990 2000

d.Soybean - Pergamino

e. Soybean - MarcosJuarez

f. Soybean - Anguil

1960 1970 1980 1990 2000

Year

20 15 0

5

10

Yield (t/ha)

20 15 10 0

0

5

10

15

Yield (t/ha)

20

25

Year

25

Year

25

Year

5

Yield (t/ha)

10 20 30

40 50 60 70

Yield (t/ha)

0

10 20 30

40 50 60 70 10 20 30 0

Yield (t/ha)

The MARS method is used to build statistical models to predict yields based on quarterly precipitation and temperature averages. The 1960–1990 years are used to train the statistical model whereas the rest of the period (1991–2005) constitutes the validation dataset. Figure 2 shows the yields from the DSSAT model and the yields predicted by the MARS models.

1960 1970 1980 1990 2000

Year

1960 1970 1980 1990 2000

Year

Fig. 2 Time series of reference yields (in tons per hectare) from the DSSAT model (“Ref.”) and predicted yields (“Pred.”) by MARS models for maize and soybean at Pergamino, Marcos Juarez and Anguil between 1960 and 2005. The dotted lines indicate the partitioning between the training dataset (1960–1990) and the validation dataset (1991–2005)

574


At first sight, the models perform well: they fit the data fairly well in the training period and the fit does not seem to drop in the validation period, apart from the case of maize at Marcos Juarez. The catastrophic yield of 1994 does not appear in the MARS model. In fact, this yield is due to low temperature values shortly after the sowing date which froze the crop. If a late sowing for 1994 had been chosen for the DSSAT model, the yield would be close to the one simulated by the MARS model. Therefore, the inaccuracy of the MARS model for this particular year is explainable and under real conditions, a farmer would certainly sow a second time or choose late sowing and obtain a yield closer to the one predicted by the MARS model. However, the MARS model also poorly fits the rest of the yield data for maize at Marcos Juarez, the trend at the end of the period being very different from what is simulated by DSSAT. Table 3 provides two statistics for the quality of the MARS models: the determination coefficient R2 is a measure of the explained variance, and the correlation coefficient is a measure of the linear dependency between the MARS predictions and the DSSAT simulations. They both clearly indicate that maize yields at Marcos Juarez are not well predicted by the MARS models in the validation period. The other predictions are found to be fairly accurate (correlations always above 0.6), especially for Pergamino. Even if they are less accurate, predictions for soybean at Anguil are robust as the summary statistics are similar for the training and the validation periods. Finally, for maize at Anguil, the correlation coefficient is above 0.7 but the determination coefficient is only 0.25, showing that the variations are quite well represented, even if the amplitudes may be over- or under-estimated (see Fig. 2c). Figure 2 shows that the MARS models present both for maize or soybean and at each site, low skill in simulating either very high or very low yields. The saturation in amplitude observed in Fig. 2 actually means that the MARS models cannot build any significant relationship based on the predictands to explain either some low- or high-peaks. This result can easily be explained by the fact that the climate–crop relationship is non-linear and that a 3-month average (especially in precipitation) may hide a large variety of high-frequency variability. Indeed, dry conditions can be observed at a critical time without being significant in a 3-month average. Therefore, this result is mainly due to seasonal climate information rather than to the statistical model and clearly indicates a limitation of the use of 3-month average climate information for crop forecasting.

Table 3 Coefficient of determination R2 and correlation ρ between the simulated yields from DSSAT ad the predicted yields by MARS models for the training period (1960–1990) and the validation period (1991–2005) at Pergamino (PE), Marcos Juarez (MJ) and Anguil (AN) for maize and soybean Crop

Site

Determination coefficient R2 Training Validation

Correlation coefficient ρ Training Validation

Maize

PE MJ AN PE MJ AN

0.58 0.68 0.45 0.83 0.80 0.58

0.76 0.83 0.67 0.91 0.89 0.67

Soybean

0.54 −0.46 0.28 0.72 0.43 0.34

0.74 −0.27 0.72 0.85 0.78 0.61


575

Table 4 Summary of MARS models to predict yields for maize and soybean at Pergamino (PE) and Anguil (AN) Crop

Site

Predictor

Range

Coefficient

Maize

PE

Soybean

AN PE

PNDJ TJFM PNDJ PDJF PJFM PNDJ PJFM

22.1◦ C >2.5 mm/day (230 mm)