# Predicting optimum crop designs using crop models and seasonal climate forecasts

D. Rodriguez, P. de Voil, D. Hudson, J. N. Brown, P. Hayman, Hélène

Marrou, H. Meinke

► To cite this version:

D. Rodriguez, P. de Voil, D. Hudson, J. N. Brown, P. Hayman, et al. Predicting optimum crop designs using crop models and seasonal climate forecasts. Scientific Reports, 2018, 8 (1), pp.1-13. 10.1038/s41598-018-20628-2. hal-01897998

HAL Id: hal-01897998

https://hal.science/hal-01897998

Submitted on 17 Oct 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

# Predicting optimum crop designs using crop models and seasonal climate forecasts

Received: 17 November 2017

Accepted: 15 January 2018

Published online: 02 February 2018

D. Rodriguez<sup>1</sup>, P. de Voil<sup>1</sup>, D. Hudson<sup>1,2</sup>, J. N. Brown<sup>3</sup>, P. Hayman<sup>4</sup>, H. Marrou<sup>5</sup> & H. Meinke<sup>6</sup>

Expected increases in food demand and the need to limit the incorporation of new lands into agriculture to curtail emissions, highlight the urgency to bridge productivity gaps, increase farmers profits and manage risks in dryland cropping. A way to bridge those gaps is to identify optimum combination of genetics (G), and agronomic managements (M) i.e. crop designs (GxM), for the prevailing and expected growing environment (E). Our understanding of crop stress physiology indicates that in hindsight, those optimum crop designs should be known, while the main problem is to predict relevant attributes of the E, at the time of sowing, so that optimum GxM combinations could be informed. Here we test our capacity to inform that "hindsight", by linking a tested crop model (APSIM) with a skillful seasonal climate forecasting system, to answer "What is the value of the skill in seasonal climate forecasting, to inform crop designs?" Results showed that the GCM POAMA-2 was reliable and skillful, and that when linked with APSIM, optimum crop designs could be informed. We conclude that reliable and skillful GCMs that are easily interfaced with crop simulation models, can be used to inform optimum crop designs, increase farmers profits and reduce risks.

Feeding the projected 9 billion people by 2050 will primarily depend on our capacity to increase food production per unit of land area, particularly in dryland cropping<fnref n="1" />, rather than incorporating new lands into agriculture²,³. Opportunities to increase yields in dryland cropping systems exist⁴,⁵, and are likely to contribute up to ca. 46% of the projected future food supply/demand gap<fnref n="1" />. Dryland cropping is typically characterized by varying levels of climate and soils variability that generate different frequencies, dynamics and intensities of water stress patterns and crop yields⁶,⁷. Amongst most broad acre dryland crops, sorghum is considered to present advantages in terms of heat and drought tolerance⁸, while recent progress in sorghum yield gains⁹ show potential for sorghum to contribute significantly to bridge expected gaps between food production and demand.

Climate variability contributes to yield gaps in three interrelated ways. (i) directly through crop physiological impacts of water and heat stresses⁷; (ii) indirectly when risk averse farmers adopt conservative management strategies that trade yield potential for reduced yield variability¹⁰; and (iii) through a 'moving target effect' where even risk neutral farmers generally opt for management options that suit an average season, but that are sub-optimal for either above or below average seasons¹¹.

In dryland sorghum cropping, growing conditions are influenced by two major factors, namely soil conditions at the time of sowing (i.e. soil fertility and soil water availability) and in-crop rainfall and its distribution during the upcoming season. While soil conditions are usually known (or at least knowable) at the beginning of the season, information about future rainfall could be obtained via skillful probabilistic seasonal climate forecasts¹². Since the early 1980s there has been a steady improvement in the skill of seasonal climate forecasts¹³. This has come in part from advances in the understanding of climate drivers and sources of predictability at a range of time scales¹⁴, but also observational and computing advances¹³. Improvements in skill have also been greatly assisted by satellite and computing technology and improvements in data assimilation and model parameterisation¹⁴. For example, a comparative assessment of accuracy and reliability of Australia's Global Circulation Model (GCM)

<footnote n="1">¹Queensland Alliance for Agriculture and Food Innovation (QAAFI), The University of Queensland, PO Box 102, Toowoomba, Queensland, 4350, Australia. ²Bureau of Meteorology, PO Box 1289, Melbourne, Victoria, 3008, Australia. ³CSIRO Agriculture and Food, 15 College Rd., Sandy Bay, Tasmania, 7005, Australia. ⁴South Australian Research and Development Institute, PO Box 397, Adelaide, South Australia, 5066, Australia. ⁵SupAgro, 2 place Viala, 34060, Montpellier, cedex 02, France. ⁶University of Tasmania, and Tasmanian Institute of Agriculture, Private Bag 98, Hobart, Tasmania, 7001, Australia. Correspondence and requests for materials should be addressed to D.R. (email: d.rodriguez@uq.edu.au)</footnote>
operational dynamic forecasting model, versus the previous statistical seasonal forecasting system, showed that the GCM based system was more reliable, consistently more accurate over a larger spatial domain, and more useful than the previous Bureau of Meteorology statistical model¹⁵.

In agriculture, seasonal climate forecasts are particularly useful when linked with crop simulation models¹⁶,¹⁷. The significance of using seasonal climate forecasts linked with crop simulation tools⁷,¹⁸,¹⁹ resides in the capacity of dynamic crop models to capture the climate - soil - crop interactions and their emerging dynamics on water and nutrient supply and demand, stress patterns and interactions determining the final yield. Recent improvements in the capacity of the APSIM (www.apsim.info) crop simulation model to model the physiology and genetics of complex adaptive traits in sorghum²⁰, provides opportunity to inform optimum combinations of genetics (G) and managements (M), known as 'crop designs' (GxM)⁷.

Both, gains in skill from the adoption of dynamic GCM forecasting systems, and improvements in crop simulation capacity²⁰,²¹, provide opportunity to develop new and more valuable climate applications. This calls for revisiting the quantification of value in the skill of existing crop modelling and climate forecasting tools to support crop design decisions and to bridge productivity gaps by providing benefits to dryland farmers world-wide. Here, the research question "What is the value of the skill in seasonal climate forecasting, to inform crop designs?" was answered by

(i) identifying a reliable and skillful seasonal climate forecast;

(ii) understanding our capacity to predict outcomes of alternative crop designs; and

(iii) using a reliable and skillful seasonal climate forecast and a crop model to quantify, in principle, the magnitude of potential changes in profits and risks from adapting crop designs to the 'where and when', i.e. location-specific and expected seasonal conditions.

## Results

**Rainfall forecasts reliability and measures of skill.** The skill of two forecast systems – the SOI phase system and the GCM POAMA-2 (see Supp. Mat.) – were compared relative to climatology. Reliability plots of three-month rainfall forecast probabilities at a lead time = 0, for above median rainfall (Fig. 1a,b), and rainfall falling in the first (Fig. 1c,d) and third (Fig. 1e,f) terciles, showed that POAMA-2 had clear advantages over the statistical SOI phase system, in terms of reduced biases, higher reliability and sharpness. The reliability for above/below median rainfall, tercile 1 and 3 rainfall were always higher for POAMA-2 i.e. in Fig. 1, *b* (slope) values were much closer to the unity, and root mean square error (RMSE) values were less than half of those for the SOI phase system. The frequency plots in the insets of Fig. 1, represent the relative fraction of grid points occurring for each forecast probability bin, and indicate a relatively sharper forecast (a greater range of probabilities issued) with POAMA-2, making it potentially more useful for decision making i.e. an improved capacity to predict extreme seasons.

Thirty-year-running percent-consistent forecasts and Brier skill score (BSS) values indicate the changing and lower skill level of the SOI phase system compared to the skill of POAMA-2 over the hindcast period 1981–2013 (Fig. 2). We also present spatial (Figs 3, 4) and seasonal (Fig. 4) variations in BSS for above/below median rainfall forecasts. In contrast to the SOI phase system, Fig. 3 shows that overall and for all the sites, the forecasts from POAMA-2 were more skillful than using climatology. For the forecasts of different seasons, POAMA-2 showed consistently higher values of BSS for above/below median rainfall forecasts across most sites compared to the SOI phase system, particularly for forecasts initialized between March and September, when POAMA-2 has its highest skill (Fig. 4). We found similar results for forecasts of rainfall falling within terciles 1 and 3 of climatology (Figures S1 and S2).

Changes in shift and dispersion between climatology and the forecasted distributions for three monthly above/below median rainfall are shown in density plots, as the relationship between the variance ratio (VR) and the absolute mean deviations (AMD) (mm) (Fig. 5 all months and locations, and Fig. S3 for each monthly forecast of three-monthly rainfall). Figure 5 shows a larger frequency of lower VR and higher AMD values for POAMA-2 than for the SOI phase system. Also, POAMA-2 consistently showed lower VR and higher AMD values for rainfall forecasts of each month of the year (Fig. S3).

Given the better performance of POAMA-2 with respect to the SOI phase system, GxExM effects and the value of the forecast skill to inform crop design was only tested using POAMA-2.

**GxExM effects on forecasts of sorghum yield.** We used POAMA-2 ensemble seasonal forecasts down-scaled to station level for the period 1981–2015 to drive multiple sets of APSIM simulations and produce simulated sorghum yields. The sets of APSIM simulations covered each possible combination of location; hybrid characteristics (i.e. maturity and tillering habit); management practices (i.e. sowing window, sowing density, row configuration, and level of nitrogen application); location (i.e. soil type); and prevailing conditions (i.e. soil water content at the time of sowing) (Table 2). The BSS and percent consistent values were calculated for sorghum yields and the distributions of BSS and percent consistent values for forecasts of above/below median, tercile 1 and tercile 3 sorghum yields are shown in Fig. 6a,b. The frequency distribution of positive (skillful) BSS values, and values of percent consistent higher than 50% varied across locations and forecasts. An analysis of variance on the values of BSS for the locations and factors in the GxExM factorial combination showed that the factor 'Location' had a highly significant (p < 0.0001) effect on the BSS (Table 3, and S1). The influence of the different factors in the GxM factorial was then tested within each location. Results showed that with the exception of soil type and sowing density at Dalby, and row configuration and nitrogen fertilization at Goondiwindi and Moree, most factors significantly affected the value of BSS (Table 3, and S1).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Reliability plots of forecast probabilities for monthly rainfall i.e. lead time = 0, for POAMA2 (a,c, and e) and the SOI phase system (b,d, and f), for above and below long-term median (a and b), and rainfall falling in the first (c and d), and third (e and f) terciles, at Moree, Goondiwindi, Dalby and Emerald Australia. In each graph, the solid line is the 1:1 relationship; b is the slope of the relationship between the observed frequency and forecasted probability; and RMSE is the root mean square error for the regression.</figcaption></figure>

We used regression trees to untangle the different GxExM effects on the value of BSS (Figs 7–10) and to answer what GxExM combinations show high or low skill values. The different G, E, and M factorial components had different importance in determining positive and high values of BSS i.e. high skill. For example, in Capella Queensland (Fig. 7) the highest values of BSS (high skill) were obtained for sorghum crops sown in October and November (Nodes 12 and 13), particularly on dry soils (Node 13). January and September sowings on high PAWC dry soils (Node 9) also showed high values of the BSS. The lowest values of BSS (no skill) were observed for December sowings (Node 2). In Dalby (Fig. 8) the highest values of BSS (high skill) were obtained for sorghum crops sown in early in September (Node 17), or with late sowings in January, of non-tillering hybrids sown at low plant densities (Node 10). In Dalby, the lowest values of BSS (no skill) were observed for December sowings (Nodes 4 and 5). In Goondiwindi (Fig. 9) the highest values of BSS (high skill) were obtained for medium and
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Thirty year running percent consistent (a to c), and Brier Skill Score (d to f) (Murphy, 1986) for the SOI phase system from 1920 to 1983 for forecasts of above and below median rainfall (a and d), rainfall falling on tercile 1 (b and e), and rainfall falling tercile 3 (c and f). Forecasts are for September to November rainfall, all sites i.e. Moree, Goondiwindi, Dalby and Emerald Australia together. The red dot is the skill for POAMA-2 for the period 1981–2013.</figcaption></figure>

<table><caption>**Table 1.** Rainfall characteristics of the locations in the study region. CV is the coefficient of variability for annual rainfall and September to February rainfall.</caption><thead><tr><td>Location</td><td>SILO station id.</td><td>Longitude</td><td>Latitude</td><td>Annual rainfall (mm)</td><td>CV</td><td>Sep-Feb rainfall (mm)</td><td>CV</td></tr></thead><tbody><tr><td>Capella, Queensland</td><td>35016</td><td>148.0</td><td>-23.1</td><td>596.1</td><td>0.34</td><td>397.1</td><td>0.40</td></tr><tr><td>Dalby, Queensland</td><td>41023</td><td>151.3</td><td>-27.2</td><td>669.8</td><td>0.25</td><td>426.2</td><td>0.32</td></tr><tr><td>Goondiwindi, Queensland</td><td>41038</td><td>150.3</td><td>-28.6</td><td>614.0</td><td>0.27</td><td>366.9</td><td>0.35</td></tr><tr><td>Moree, New South Wales</td><td>53027</td><td>149.9</td><td>-29.5</td><td>584.3</td><td>0.28</td><td>339.1</td><td>0.38</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Brier Skill Score for above and below median, tercile 1 and tercile 3 three monthly rainfall (lead time = 0) from January to December, using POAMA-2 and the SOI phase system, at all three locations (All), Capella (Cap), Dalby (Dalb), Goondiwindi (Gnd), and Moree (Mre) Australia.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Brier Skill Scores for forecasts of above/below median three-monthly rainfall forecasted using POAMA2 and the SOI phase system for: all three locations (All), and for Capella (Cap), Dalby (Dalb), Goondiwindi (Gnd), and Moree (Mre) Australia. The month indicated on the figure corresponds to the first month of the forecast (i.e. for "Oct" the forecast is initialized on the 1st Oct and is assessed for the OND season).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Relationship between the variance ratio i.e. a measure of dispersion, and the absolute mean deviation (mm) i.e. a measure of shift, for forecasts of above median three-monthly rainfall using POAMA2 and the SOI phase system. The figure is for all months and locations in the study i.e. Capella, Dalby, and Goondiwindi, and Moree Australia.</figcaption></figure>
<table><caption>**Table 2.** Location and soil denominator, soil type, and soil water characteristics at the locations in the study region. Plant available water (PAW) is the difference between the drainage upper limit and the crop lower limit; Plant available water capacity (PAWC) is the difference between the drainage upper limit and wilting point.</caption><thead><tr><td>Location</td><td>Soil type</td><td>PAW (mm)</td><td>PAWC (mm)</td></tr></thead><tbody><tr><td>Capella, high</td><td>Black clay Vertosol</td><td>287</td><td>415</td></tr><tr><td>Capella, medium</td><td>Grey medium clay Vertosol</td><td>189</td><td>274</td></tr><tr><td>Capella, low</td><td>Black light clay Vertosol</td><td>146</td><td>146</td></tr><tr><td>Dalby, high</td><td>Black clay Vertosol</td><td>308</td><td>400</td></tr><tr><td>Dalby, medium</td><td>Black medium clay Vertosol</td><td>253</td><td>344</td></tr><tr><td>Dalby, low</td><td>Grey light clay Vertosol</td><td>175</td><td>274</td></tr><tr><td>Goondiwindi, high</td><td>Grey clay Vertosol</td><td>283</td><td>313</td></tr><tr><td>Goondiwindi, medium</td><td>Grey medium clay Vertosol</td><td>214</td><td>361</td></tr><tr><td>Goondiwindi, low</td><td>Grey light clay Vertosol</td><td>179</td><td>310</td></tr><tr><td>Moree, high</td><td>Black clay Vertosol</td><td>299</td><td>254</td></tr><tr><td>Moree, medium</td><td>Grey medium clay Vertosol</td><td>220</td><td>304</td></tr><tr><td>Moree, low</td><td>Black light clay Vertosol</td><td>106</td><td>124</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Brier Skill Scores (a) and Percent consistent values (b) for forecasts of above or below median, tercile 1, and tercile 3 sorghum yields simulated using a factorial combination of genetic, management and site conditions, using POAMA- at Capella, Dalby, and Goondiwindi, and Moree Australia.</figcaption></figure>

late sorghum crops sown in October and November (Node 11) or early hybrids sown in November (Node 10). In Goondiwindi, the lowest values of BSS (no skill) were observed for December sowings (Nodes 3). In Moree (Fig. 10) the highest values of BSS (high skill) were obtained for sowings in October and high sowing densities (Node 13), or low sow densities of high tillering hybrids (Node 12). In Moree, the lowest values of BSS (no skill) were observed for December sowings particularly on dry soils (Nodes 3), or October sowings of no tillering late hybrids.

In general, the results from the regression trees across all locations followed a similar pattern: higher values of BSS (high skill) were observed when the GxExM factor combination resulted in crop designs that were highly dependent on seasonal conditions, i.e. in-crop rainfall.

**Value of APSIM and POAMA-2 to inform crop designs (GxM).** The value of linking APSIM and a GCM (POAMA-2) to optimize crop designs for each individual year in the hindcast series requires time series of simulated sorghum yields for three distinct different crop design (GxM) strategies. Reflecting current farmer management (GxM<sub>f</sub>) we used a single and static hybrid by management combination optimized for most years (GxM<sub>optS</sub>); and a dynamic GxM choice optimized for each individual year in response to the seasonal climate forecast (GxM<sub>optSCF</sub>). The farmer GxM<sub>f</sub> and optimum GxM<sub>optS</sub> managements do not change during the hindcast period 1981 to 2015. In contrast, GxM<sub>optSCF</sub> is adjusted each year in response to the climate forecast, equations (1) and (2). At each location, the three-time series of simulated yields were converted to an economic time series of profit ($/ha) by multiplying the yield (t/ha) by a price ($/t) to calculate annual income ($/ha) and subtracting a
<table><caption>**Table 3.** Levels of significance from ANOVA for variables related to genotype characteristics (G), management practices (M) and location (E), for all the locations together (All sites); and Capella, Dalby, Goondiwindi, and Moree individually, for the GxMxE combinations showing Brier Skill Scores in the upper 75% of the values. Significance codes: ***0; **0.001; **0.01; and *0.05; ns = not significant. na = not applicable.</caption><thead><tr><th></th><th>All sites</th><th>Capella</th><th>Dalby</th><th>Goondiwindi</th><th>Moore</th></tr></thead><tbody><tr><td>Location</td><td>***</td><td>na</td><td>na</td><td>na</td><td>na</td></tr><tr><td>Sowing window</td><td>***</td><td>***</td><td>***</td><td>***</td><td>***</td></tr><tr><td>Soil type</td><td>***</td><td>***</td><td>ns</td><td>***</td><td>***</td></tr><tr><td>Plant available water</td><td>***</td><td>***</td><td>***</td><td>***</td><td>***</td></tr><tr><td>Hybrid maturity</td><td>***</td><td>***</td><td>***</td><td>***</td><td>***</td></tr><tr><td>Tillering type</td><td>***</td><td>***</td><td>**</td><td>*</td><td>*</td></tr><tr><td>Density</td><td>***</td><td>***</td><td>ns</td><td>***</td><td>***</td></tr><tr><td>Row configuration</td><td>***</td><td>**</td><td>***</td><td>ns</td><td>ns</td></tr><tr><td>Nitrogen fertilization</td><td>***</td><td>***</td><td>***</td><td>ns</td><td>ns</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Regression tree on the Brier Skill Scores for forecasts of above median sorghum yields simulated using a factorial combination of genetic, management and site conditions, using POAMA-2 at Capella Australia. References: Window indicates the sowing month; Soil refers to the soil high (h), medium (m) or low (l) plant available water capacity; ISW refers to the soil plant available water content (v v⁻¹) at the time of sowing; and N to kg N ha⁻¹ applied at sowing. In the boxplots n refers to the number of node simulations and the red dashed line indicates a Brier Skill Score value of 0.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Regression tree on the Brier Skill Scores for forecasts of above or below median sorghum yields simulated using a factorial combination of genetic, management and site conditions, using POAMA2 at Dalby Australia. References: Window indicates the sowing month; Tillers refers to the tillering type of the hybrid where 0 is no tillering and 4 is high tillering; ISW refers to the soil plant available water content (v v⁻¹) at the time of sowing; Density to the number of plants sown m⁻²; Cultivar is the maturity type of the hybrid, i.e. early, medium, or late; and N to kg N ha⁻¹ applied at sowing. In the boxplots n refers to the number of node simulations and the red dashed line indicates a Brier Skill Score value of 0.</figcaption></figure>
<table><caption>**Table 4.** Mean profits from a simulation using farmers current practice and climatology (Farmers), perfect knowledge (PK) i.e. optimized crop designs using climatology, optimized crop design using POAMA-2 (Optimised); and value i.e. mean changes in profit and down side risk i.e. likelihood of a profit lower than $600 ha⁻¹ for, farmers current practice (Value<sub>opt</sub>). Optimized crop design using climatology 1981–2015 and (Value<sub>optS</sub>) and a POAMA-2 seasonal climate forecast (Value<sub>optSCF</sub>). Value<sub>optS</sub> = difference in profit ($ ha⁻¹) and down side risk (DSR), between simulation of current farmers’ hybrid by management combination (Value<sub>opt</sub>) and a static (every year the same) optimized hybrid by management combination. Value<sub>optSCF</sub> = difference in profit ($ ha⁻¹), between Value<sub>optS</sub> and the dynamically (every year different) optimized hybrid by management combination informed by the POAMA-2 seasonal climate forecasts. Value<sub>PK</sub> = Value of having perfect knowledge i.e. optimum crop design using observed climatology.</caption><thead><tr><th rowspan="3">Soil type (PAWC)</th><th colspan="6">Mean value to inform GxM with respect to a static management</th><th colspan="3">DSR (%)</th></tr><tr><th colspan="4">Profit ($ ha<sup>-1</sup>)</th><th rowspan="2">Value<sub>optSCF</sub></th><th rowspan="2">Value<sub>PK</sub></th><th rowspan="2">Farmers</th><th rowspan="2">Optimized</th><th rowspan="2">Value<sub>optS</sub></th></tr><tr><th>Farmers</th><th>Optimized</th><th>Value<sub>optS</sub></th><th>Value<sub>optSCF</sub></th></tr></thead><tbody><tr><td rowspan="3">Capella</td><td>High</td><td>1108</td><td>1260</td><td>152</td><td>3</td><td>335</td><td>0</td><td>0</td><td>No risk</td></tr><tr><td>Medium</td><td>748</td><td>824</td><td>77</td><td>3</td><td>224</td><td>12</td><td>0</td><td>-100</td></tr><tr><td>Low</td><td>544</td><td>600</td><td>56</td><td>4</td><td>219</td><td>68</td><td>62</td><td>-9</td></tr><tr><td rowspan="3">Dalby</td><td>High</td><td>1127</td><td>1337</td><td>210</td><td>13</td><td>353</td><td>0</td><td>0</td><td>No risk</td></tr><tr><td>Medium</td><td>1048</td><td>1241</td><td>194</td><td>17</td><td>351</td><td>0</td><td>0</td><td>No risk</td></tr><tr><td>Low</td><td>795</td><td>913</td><td>118</td><td>12</td><td>288</td><td>6</td><td>3</td><td>-50</td></tr><tr><td rowspan="3">Goondiwindi</td><td>High</td><td>866</td><td>1092</td><td>226</td><td>16</td><td>432</td><td>6</td><td>0</td><td>-100</td></tr><tr><td>Medium</td><td>841</td><td>1011</td><td>170</td><td>63</td><td>345</td><td>3</td><td>0</td><td>-100</td></tr><tr><td>Low</td><td>678</td><td>793</td><td>115</td><td>6</td><td>312</td><td>34</td><td>6</td><td>-82</td></tr><tr><td rowspan="3">Moree</td><td>High</td><td>1025</td><td>1226</td><td>202</td><td>23</td><td>406</td><td>0</td><td>0</td><td>No risk</td></tr><tr><td>Medium</td><td>814</td><td>962</td><td>148</td><td>32</td><td>370</td><td>9</td><td>0</td><td>-100</td></tr><tr><td>Low</td><td>373</td><td>427</td><td>54</td><td>19</td><td>210</td><td>89</td><td>86</td><td>-3.3</td></tr></tbody></table>

profits were 21% in also 26% of the years (Table S2). When we calculated the value of skill in seasonal climate forecasting relative to an improved static management involving higher levels of investment in nitrogen fertilizers, the additional value of the new climate information was smaller, on average 17 AU$ ha⁻¹, i.e. Value<sub>optSCF</sub> in Table 4. In Table 4 we also show the value of perfect knowledge (Value<sub>PK</sub>), this is, the difference between an optimized crop design and the static farmers’ management, calculated using observed climatology. From comparing Value<sub>PK</sub>, and Value<sub>optS</sub> (Table 4) we can conclude that present value in the skill of seasonal climate forecasting falls approximately mid-way between ‘no skill’ and perfect knowledge of the future climate.

Discussion

Inherent to dryland cropping is a high level of season-to-season and within-season climate variability¹⁴,²². This, together with different G and M combinations results in highly contrasting stress environments and crop yields⁶,⁷. Australia has a long track record of valuable developments in climate sciences and applications²³, such as the development and use of seasonal climate forecasts to inform likely seasonal conditions and practice change¹⁵,²³, see also Climate Kelpie at http://www.climatekelpie.com.au. However, adoption remains low due to, (i) the perceived low value of the existing skill in the information of seasonal climate forecasts²⁴; (ii) the complexities associated with the multiple interactions between factors when managing biological systems (i.e. climate, soil and crop interactions, and their effect on the skill and value of crop yield forecasts)²³; and (iii) the challenge of understanding and communicating probabilistic information²⁴, especially by risk averse farm managers and consultants¹⁰.

Here we used measures of forecast reliability, skill, shift and dispersion to identify a reliable and skillful seasonal climate forecast for rainfall. We then used the most reliable and skillful forecast in combination with a validated crop simulation model to capture the multiple interactions between sorghum hybrid, management and environmental conditions. Finally, we calculated the value of skill in seasonal climate forecasting to inform optimum crop designs for sorghum in Australia.

Results showed that the forecasts produced by Australia's operational forecasting system POAMA-2²⁵ proved to be more reliable, sharper and consistently more accurate than the SOI phase system²⁶. The extent of the difference between both systems was substantial; we know of no previous attempts to compare them. Results show that the dynamic climate model outperformed the statistical climate model across all tested indices (Figs 1 and 2). In general, the SOI phase system was worse than climatology at predicting rainfall. This is likely a consequence of progressive skill degradation over time since the system was first developed in the mid 1990s²⁶. Being a static, statistical system based on historical relations between the SOI and future rainfall, it is likely to have suffered from degradation due to climate change and low frequency climate fluctuations that have affected the historical climate record since²⁷. Significant improvements in sharpness, reliability and accuracy from the previous Bureau of Meteorology (BOM) statistical models²⁸ to the current BOM's dynamic models²⁵ were reported before¹⁵. Once challenges such as data initialization and resolution are addressed, dynamical models are logically able to outperform
their simpler statistical counterparts for a number of reasons. Firstly, dynamic models can incorporate a number of climate drivers beyond ENSO such as the Indian Ocean Dipole and monsoon variability and their complex interrelationships. Secondly, climate change and its effect on climate drivers are incorporated into dynamical models whereas statistical approaches are limited to analogues from the past²⁹. Also, decadal variations in climate interactions require a long historical record for differing phases to be accounted for³⁰, particularly the variations in predictive skill across decades²⁷. BSS values for above/below median rainfall, and rainfall falling within tercile 1 or 3, varied across sites (Fig. 3) and within sites during the year (Fig. 4). This is consistent with observations of spatial and temporal variations in the influence of ENSO across eastern Australia¹⁴. An important observation was the significantly higher level of skill in POAMA-2 between March and September across all tested locations. Although the locations in this study are mostly summer rainfall dominant, this is at a time of the year when some of the most profitable crops are grown in the region e.g. wheat and chickpeas.

When we calculated BSS values for a simulated factorial combination of GxM factors on sorghum yields, we observed a large variation in BSS and percent consistent values for forecasts of above/below median, tercile 1 and tercile 3 sorghum yields (Fig. 6a,b). The large variability in BSS values from the GxMxE factorial combination was expected, as the influence of seasonal climate forecasts on sorghum yield is likely to vary with the reliance of the crop design on in-crop rainfall and seasonal conditions³¹. This is consistent with others³², who found that sampling variances increase with increasing forecast accuracy, and with decreasing climatological event probability. Our results showed that location, site, hybrid and management practices all significantly affected the value of BSS, and that as indicated above, the highest values of BSS were observed for crop designs that were highly reliant on in-crop rainfall. The capacity of soils to buffer impacts climate variability and to modify the skill and value of a forecast was shown before³¹. The influence of site, hybrid, and management combinations on the economic value of using a crop model and GCM forecasts is discussed below.

Reliability and skill in a seasonal climate forecast is essential, although it does not guarantee that the information has value, or is of sufficient importance to change the behavior of mostly risk averse decision makers such as dryland farmers¹⁰,³³. Previous analyses showed that up to 20% increases in profit and or up to 35% reductions in risk were possible when tactical adjustments to crop nitrogen management were informed using the SOI phase system¹⁷. More recently, in silico estimates of benefits from environment specific crop designs e.g. associated with increasing plant densities and sowing low tillering sorghum hybrids, ranged between 0.22 and 0.41 t ha⁻¹ ca. 56–104 AU$ ha⁻¹ at average sorghum prices (254 AU$ t⁻¹), together with slight reductions in down side risk⁷. Here we showed that the value, in terms of increases in profits and reductions in downside risks, from identifying optimum GxM combinations using a crop simulation model and a seasonal climate forecast was significant across all tested locations and soil conditions i.e. soil PAWC values (Table 4). The value in skill depended on the baseline for the comparison: When current farmers’ practice was used as the baseline, linking APSIM sorghum and POAMA-2 increased average profits by 143 AU$ ha⁻¹ and reduced or even eliminated down side risk. When the baseline for the comparison was the highest yielding, static hybrid-by-management combination (Table 4 and S2), the actual value of the additional climate information was on average 17 AU$ ha⁻¹, which compares to the benefits derived from Australia’s sorghum breeding over the last thirty years⁹ i.e. 2.1% per year, or 44 kg ha⁻¹ year⁻¹. These results indicate that even though the value of the additional climate information might seem small (Valueₐ₉₄₅), its magnitude compares well with that derived from much larger and better funded breeding programs. Much larger benefits (Valueₐ₉₅₆) might be realized when using such insights in discussions with farmers on benefits and risk from increasing investments in dryland cropping¹⁰ to sustainably bridge productivity and profit gaps in dryland cropping i.e. the difference between Profitₐ₉₇₇ and Profitₐ₉₇₈ or Profitₐ₉₇₈ (Table 4).

Given the chaotic nature of the atmosphere it is impossible to know exactly how it will evolve beyond a few days³⁴, requiring the use of probabilistic forecasts and analyses of changes in probability distribution functions, with respect to climatology¹²,¹⁸,³⁴. However, given the hypothetical case that we would have perfect knowledge about the future, we also estimated that the present value in seasonal climate forecasting falls mid-way from having perfect knowledge (Valueₚₖ in Table 4). Clearly, communicating probabilistic information requires representing the complete distribution of likely outcomes from a change in practice against the counterfactual of ‘no action’, and accepting that pay-offs will be realized long term, rather than every single occasion the forecast is used. This is important, given that the wide range of individual circumstances affects farmers’ levels of risk aversion and investment capacity. The results presented here indicate that improved agronomic practice and increased investments in fertilizers can lift average profits by up to 82%, while climate forecasts information can contribute up to 12%, a value that is on par with plant breeding.

## Conclusions

We conclude that reliable and skillful dynamic GCM models, interfaced with validated crop simulation models, can now be used to inform optimum crop designs to increase farmers’ profits and reduce risks. Australia’s Bureau of Meteorology current POAMA-2 model outperformed the statistical SOI phase system across all tested measures of reliability, skill, shift and dispersion, particularly during autumn and winter. The expected release of BOM’s new higher resolution and more sophisticated ACCESS-S1 seasonal climate forecast system early during 2018 is likely to increase further the value of climate information when linked with crop simulation models like APSIM. To achieve those gains, improvements in downscaling techniques and real-time access to outputs from BOM’s seasonal climate forecasts will be required. Further, efforts need to continue that makes such information available to decision makers in a form that is understandable and useable.

## Methods

**Climate data and seasonal climate forecasting tools.** Long term climate records, the SOI phase system²⁶ and a hindcast from the operational Australian seasonal climate forecasting system (POAMA-2)²⁵ were assembled for four sites in Australia’s northern grains region i.e. Capella, Dalby and Goondiwindi in Queensland,
and Moree in New South Wales (Table 1). Weather data included both station data from SILO (1889–2016) (https://www.longpaddock.qld.gov.au/silo/), and gridded (2.5 degree) data (1901–2013) from GPCC (https://www.esrl.noaa.gov/psd/data/gridded/data.gpcc.html). At each location, “0-lead” rainfall forecasts for the next ninety days were made at the beginning of each month using both POAMA-2 and the SOI phase system for the common hindcast period of 1981–2015.

The SOI phase system²⁶ was used to develop forecast distributions from analogue years of observed data, a cross validated probability of falling into each category was calculated as in¹⁷. The SOI phase system uses five fixed phases derived from consecutive monthly values of the SOI normalized with respect to the base period 1887–1989. The five phases of the SOI are rapidly rising, consistently positive, neutral or near zero, rapidly falling and consistently negative. Analog years corresponding to each of the SOI phases were obtained from the Long Paddock website (www.longpaddock.qld.gov.au). The analogue sets used to create the SOI phase distributions were taken from the period 1901–2013.

POAMA-2 is a global ensemble seasonal forecast system, comprising a coupled ocean-atmosphere model and data assimilation systems for the initialization of the ocean, land and atmosphere²⁵. The forecasts are initialized with observational data available at the start of the forecast (to provide the best representation of the current state of the climate system) and are then run forward in forecast mode for nine months. Multiple forecasts are run for a given start time, called an ensemble, in order to provide an indication of the uncertainty of the future evolution of the climate system i.e. the likelihood of future conditions is presented as a forecast distribution. The ensemble for a given start time consists of thirty-three forecasts that differ only in their initial conditions. POAMA-2 forecasts are made every 5 days and in order to expand the number of ensemble members, the outputs from three consecutive runs were used to create a 99-member ensemble (this approach is also taken by the Bureau of Meteorology for their climate outlooks). POAMA-2 forecasts were made for four locations each falling into a separate 2.5-degree grid cell. To assess the value of using POAMA-2 seasonal forecasts, we use a set of forecasts that have been run retrospectively over a period in the past, called a hindcast set. The hindcast period is 1981–2015 and the 99-member ensemble forecasts issued on the first day of every month in that period were used in this study. POAMA-2 forecasts were assembled in two forms, (i) a 3-month total rainfall forecast made on the first day of every month; and (ii) a downscaled daily time series of rainfall, maximum and minimum temperatures and solar radiation, suitable to be used with a crop simulation model³⁵.

**APSIM modelling.** The open source code validated APSIM sorghum model²⁰ was used to simulate a comprehensive factorial combination of hybrid characteristics (G), management factors and site-soil conditions across four locations in Queensland and New South Wales, Australia. APSIM sorghum was developed and parameterized using local field experimentation and cultivars, details on model equations and validation results are provided with the supplementary material, and can be also accessed at http://www.apsim.info/APSIM Validation/ Main.aspx

Simulations for the period 1981–2015 using observed station data (n = 1) and the downscaled daily time series forecast (n = 99) from POAMA were conducted for monthly sowing windows between September and January; on three different soil types per location (Table 2); and four values of initial soil water (ISW) at the time of sowing (0.2, 0.4, 0.6, and 0.8 v/v). Agronomic management factors included two row configurations (solid 1 m rows, and single skip)³⁶, four levels of plant density (3.5, 5, 6.5, 8 pl m⁻²), and three levels of nitrogen fertilization (50, 100 and 150 kg N ha⁻¹). Genotype characteristics included three levels of maturity (early, medium and late), and three levels of tillering (no tillering, medium and high tillering) as in⁷. The model was run over the 1981–2015 (34 years) hindcast series, resulting in a total number of GxExM combinations of 176,256,000 crop year simulations.

**Measures of skill.** The skill of both forecasting systems was calculated on both 3-month gridded rainfall (GPCC) and simulated sorghum yields. From the rainfall forecast distributions, probabilities of above median, tercile 1 and 3 rainfall, and a binary outcome i.e. correct or incorrect, for each forecast were calculated. Measures of skill included the reliability diagram³⁷, the Brier Skill Score (BSS, equations (S3 and S4))³⁸, percent consistent rates, and measures of shift in the mean i.e. absolute mean deviation (AMD, equation (S5)), and dispersion i.e. the variance ratio (VR, equation (S6)) of the sample variances¹⁸,³⁹.

Percent consistent compares how often the forecast favored a particular outcome and how often that outcome was realized¹⁵. Reliability diagrams indicate errors associated with the issued probabilities. In the reliability plots the y-axis is the relative observed frequency and the x-axis is the forecast probability. The solid line in the reliability plots indicates perfect reliability. The reliability diagram is accompanied by a histogram which indicates the sample size in each forecast probability bin. The histogram can be used to indicate possible sampling issues associated with each bin (e.g. too few samples for a robust result), as well as the sharpness of the forecast system. If all the forecasts fell in the model climatology probability bin, then the system would have no sharpness (sharpness is the tendency to forecast extreme values). Sharp forecasts are useful for decision making, assuming they are reliable. The Brier Skill Score indicates the proportional improvement of the probabilistic forecasts from a given system over using climatological forecasts. All statistical analyses were calculated using the R software³⁸.

**Value of linking APSIM and POAMA-2 to inform crop designs (GxM).** The value of using POAMA-2 to inform more profitable and resilient crop designs was derived for the subset of GxExM combinations producing positive BSS values for forecasts of above/below median sorghum yields. A search algorithm was developed that included the following steps (i) a matrix of average profits and down side risk for best farmers' practice and for a factorial combination of GxM (hybrid maturity and tillering type, plant density, row configuration and fertilization rate) was created; (ii) GxM combinations having values of BSS higher than zero (skill) for predicted profits were kept; (iii) For each year in the hindcast the GxM alternative showing the highest profit were then
selected; (iv) Value was then calculated in terms of changes in average profits and down side risk i.e. likelihood of a profit lower than 600 AU$ ha⁻¹, between best farmers' practice, i.e. a static though locally recommended best agronomic management, and the selected GxM combination predicted using the POAMA-2 forecast for each year in the hindcast series, 1981–2015 (Valueₐₗₕ); and relative to the simulated most frequent optimum static GxM strategy (ValueₐₗₕC̈F̈) (Table S2). In this analysis Valueₐₗₕ represents both gains from improved crop design and the use of climate information, while ValueₐₗₕC̈F̈ represents the actual value of the new climate information. As a reference point, value was also calculated for the hypothetical situation that the future climate was known i.e. perfect knowledge (Valueₚₖ). Valueₚₖ was calculated as the difference in profit between the optimized crop design and farmers' practice using observed climatology. To account for opportunity and fixed costs, a threshold value of 600 AU$ ha⁻¹ to calculate down side risk was chosen after consultation with farmers across the locations in the study. Profits were calculated for median sorghum prices over the last 10 years (2007 to 2016) i.e. 254 AU$ t⁻¹. Variable costs included insurance i.e. 1% of the gross income, while the fertilizer cost was 30 cents kg⁻¹ urea (http://agmargins.net.au). Variable costs for failed crops were then 161 AU$ ha⁻¹ (excluding harvest costs), and 211 AU$ ha⁻¹ for harvested crops.

**Data and model availability.** APSIM is an open source cropping systems model available at www.apsim.info.

**References**

1. Keating, B. A., Herrero, M., Carberry, P. S., Gardner, J. & Cole, M. B. Food wedges. Framing the global food demand and supply challenge towards 2050. *Global Food Security* 1–8 https://doi.org/10.1016/j.gfs.2014.08.004 (2014).

2. Godfray, H. C. J. & Garnett, T. Food security and sustainable intensification. *Philosophical Transactions of the Royal Society B: Biological Sciences* **369**, 20120273–20120273 (2014).

3. Garnett, T. et al. Sustainable Intensification in Agriculture: Premises and Policies. *Science* **341**, 33–34 (2013).

4. Lobell, D. B., Cassman, K. G. & Field, C. B. Crop Yield Gaps: Their Importance, Magnitudes, and Causes. *Annu. Rev. Environ. Resour.* **34**, 179–204 (2009).

5. van Ittersum, M. K. et al. Yield gap analysis with local to global relevance - A review. *Field Crops Research* **143**, 4–17 (2013).

6. Chauhan, Y. S., Solomon, K. F. & Rodriguez, D. Characterization of north-eastern Australian environments using APSIM for increasing dryland maize production. *Field Crops Research* **144**, 245–255 (2013).

7. Hammer, G. L. et al. Crop design for specific adaptation in variable dryland production environments. *Crop Pasture Sci.* **65**, 614–624 (2014).

8. Borrell, A. K. et al. Drought adaptation of stay-green sorghum is associated with canopy development, leaf anatomy, root growth, and water uptake. *Journal of Experimental Botany* **65**, 6251–6263 (2014).

9. Potgieter, A. B. et al. Yield trends under varying environmental conditions for sorghum and wheat across Australia. *Agricultural and Forest Meteorology* **228–229**, 276–285 (2016).

10. Monjardino, M., McBeath, T. M., Brennan, L. & Llewellyn, R. S. Are farmers in low-rainfall cropping regions under-fertilising with nitrogen? A risk analysis. *Agricultural Systems* **116**, 37–51 (2013).

11. Meza, F. J., Hansen, J. W. & Osgood, D. Economic Value of Seasonal Climate Forecasts for Agriculture: Review of Ex-Ante Assessments and Recommendations for FutureResearch. *J. Appl. Meteor. Climatol.* **47**, 1269–1286 (2008).

12. Hansen, J., Baethgen, W., Osgood, D., Ceccato, P. & Ngugi, R. K. Innovations in Climate Risk Management: Protecting and Building Rural Livelihoods in a Variable and Changing Climate. *Journal of Semi-Arid Tropical Agricultural Research* **4**, 1–38 (2007).

13. Bauer, P., Thorpe, A. & Brunet, G. The quiet revolution of numerical weather prediction. *Nature* **525**, 47–55 (2015).

14. Risbey, J. S., Pook, M. J., McIntosh, P. C., Wheeler, M. C. & Hendon, H. H. On the Remote Drivers of Rainfall Variability in Australia. *Mon. Wea. Rev.* **137**, 3233–3253 (2009).

15. Charles, A. N., Duell, R. E. & Wang, X. *Seasonal Forecasting for Australia using a Dynamical Model: Improvements in Forecast Skill over the Operational Statistical Model. Australian Meteorological and Oceanographic Journal* **65**, 3–4 (AUSTRALIAN ... , 2015).

16. Hansen, J. Linking dynamic seasonal climate forecasts with crop simulation for maize yield prediction in semi-arid Kenya. *Agricultural and Forest Meteorology* **125**, 143–157 (2004).

17. Hammer, G. L., Holzworth, D. & Stone, R. The value of skill in seasonal climate forecasting to wheat crop management in a region with high climatic variability. *Aust J Agric Res* **47**, 717–737 (1996).

18. Potgieter, A. B., Everingham, Y. L. & Hammer, G. L. On measuring quality of a probabilistic commodity forecast for a system that incorporates seasonal climate forecasts. *Int. J. Climatol.* **23**, 1195–1210 (2003).

19. Crane, T. A. et al. Forecast Skill and Farmers' Skills: Seasonal Climate Forecasts and Agricultural Risk Management in the Southeastern United States. *Wea. Climate Soc.* **2**, 44–59 (2010).

20. Hammer, G. L. et al. Adapting APSIM to model the physiology and genetics of complex adaptive traits in field crops. *Journal of Experimental Botany* **61**, 2185–2202 (2010).

21. Holzworth, D. P. et al. APSIM - Evolution towards a new generation of agricultural systems simulation. *Environmental Modelling & Software* **62**, 327–350 (2014).

22. van Dijk, A. I. J. M. et al. The Millennium Drought in southeast Australia (2001–2009): Natural and human causes and implications for water resources, ecosystems, economy, and society. *Water Resour. Res.* **49**, 1040–1057 (2013).

23. Hammer, G. L. et al. Advances in application of climate prediction in agriculture. *Agricultural Systems* **70**, 515–553 (2001).

24. Hayman, P., Crean, J., Mullen, J. & Parton, K. How do probabilistic seasonal climate forecasts compare with other innovations that Australian farmers are encouraged to adopt? *Aust. J. Agric. Res.* **58**, 975–10 (2007).

25. Hudson, D., Marshall, A. G., Yin, Y., Alves, O. & Hendon, H. H. Improving Intraseasonal Prediction with a New Ensemble Generation Strategy. *Mon. Wea. Rev.* **141**, 4429–4449 (2013).

26. Stone, R., Hammer, G. L. & Marcusen, T. Prediction of global rainfall probabilities using phases of the Southern Oscillation Index. *Nature* **384**, 252–255 (1996).

27. Zhao, M., Hendon, H. H., Alves, O., Liu, G. & Wang, G. Weakened Eastern Pacific El Niño Predictability in the Early Twenty-First Century. *J. Climate* **29**, 6805–6822 (2016).

28. Drosdowsky, W. & Chambers, L. Near-global sea surface temperature anomalies as predictors of Australian seaosonal rainfall. *J. Climate* **14**, 1677–1687 (2001).

29. Power, S. B. & Kociuba, G. The impact of global warming on the Southern Oscillation Index. *Clim Dyn* **37**, 1745–1754 (2010).

30. Kirtman, B. P. & Schopf, P. S. Decadal variability in ENSO predictability and prediction. *J. Climate* **11**, 2804–2822 (1998).

31. Wang, E., Cresswell, H., Xu, J. & Jiang, Q. Capacity of soils to buffer impact of climate variability and value of seasonal forecasts. *Agricultural and Forest Meteorology* **149**, 38–50 (2009).

32. Wilks, D. S. Sampling distributions of the Brier score and Brier skill score under serial dependence. *Quarterly Journal of the Royal Meteorological Society* **136**, 2109–2118 (2010).
33. Meinke, H. & Stone, R. Seasonal and inter-annual climate forecasting: The new tool for increasing preparedness to climate variability and change in agricultural planning operations. *Climatic Change* **70**, 221–253 (2005).

34. Goddard, L. et al. Current approaches to seasonal to interannual climate predictions. *Int. J. Climatol.* **21**, 1111–1152 (2001).

35. McIntosh, P. C. & Brown, J. Calibration and bias correction of seasonal climate forecasts for use in agricultural models. 1–28, https://doi.org/10.4225/08/5910c03c4d437 (2017).

36. Whish, J. et al. Modelling the effects of row configuration on sorghum yield reliability in north-eastern Australia. *Aust. J. Agric. Res.* **56**, 11–13 (2005).

37. Bröcker, J. & Smith, L. A. Increasing the Reliability of Reliability Diagrams. *Wea. Forecasting* **22**, 651–661 (2007).

38. Murphy, A. H. A new decomposition of the Brier score: Formulation and interpretation. *Mon. Wea. Rev.* **114**, 2671–2673 (1986).

39. Fawcett, F. & Stone, R. A comparison of two seasonal rainfall forecasting systems for Australia. *Australian Meteorological and Oceanographic Journal* **60**, 15–24 (2010).

**Acknowledgements**

This work was supported by the Grains Research and Development Corporation through the project “High yielding agronomy for wheat and sorghum in the northern grains region” UQ00074. PTH time was funded by the Australian Government through the RR&D for Profit project “Improved use of SCF to increase farmer profitability”. We thank Prof Graeme Hammer for his comments on an early version of this manuscript, and the reviewers who contributed to improve the quality of the manuscript.

**Author Contributions**

D.R. wrote the manuscript; P.de.V. run APSIM and POAMA-2, prepared the figures and contributed to the writing of the manuscript; D.H., J.B., P.H., H.Ma. and H.Me. contributed to the writing of the manuscript.

**Additional Information**

**Supplementary information** accompanies this paper at https://doi.org/10.1038/s41598-018-20628-2.

**Competing Interests:** The authors declare that they have no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_11.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018