Probabilistic modeling of crop-yield loss risk under drought: a spatial showcase for sub-Saharan Africa

We investigated three staple crops (i.e. maize, millet, and sorghum), widely grown and adapted to the environmental conditions of SSA, for their specific response to drought events. The crop yields were simulated throughout SSA during 1980–2012 period, at 0.5° spatial resolution using the environmental policy integrated climate (EPIC) model version extended with a calibration module i.e. EPIC⁺ [44]. EPIC simulates crop growth processes in daily time steps using weather, soil, land use, and crop management parameters [45]. EPIC⁺ develops a spatially explicit version of original EPIC in a Python framework. This allows the extension of the EPIC application to larger scales by dividing the region of study into grids based on a specified resolution (here 0.5°) and executing EPIC on each grid cell.

Additionally, EPIC⁺ uses the sequential uncertainty fitting (SUFI-2) algorithm [46] for automatic calibration and uncertainty assessment of the model. According to SUFI-2 algorithm, uncertainty is quantified as the 95% prediction uncertainty (95PPU), calculated at the 2.5% and 97.5% levels of the variable's cumulative distribution using Latin hypercube sampling in the parameters' space. We used two statistics—the p-factor and the r-factor—to quantify the goodness of model calibration and the uncertainty level [47]. The p-factor, the fraction of the observed data bracketed by the 95PPU band, varies from 0 to 1, where 1 indicates 100% bracketing of the observed data within the model PPU. Values around 0.5 are often considered acceptable for crop simulations [48]. In contrast, the r-factor is the ratio of the average width of the 95PPU band to the standard deviation of the measured variable. The ideal value for the r-factor is 0, with an acceptable practical value of around 2 for crop yields [4]. The SUFI-2 algorithm endeavors to achieve a high p-factor while keeping the r-factor as small as possible [48].

The simulated maize yields are obtained from a study conducted by Kamali et al [20] in which the resulted gridded yields were aggregated to country level for model calibration. Here, we implemented the same procedure for simulation runs and model calibration for the case of sorghum and millet crops across SSA. The site-specific input data included longitude, latitude, slope, elevation (DEM), climate, soil, crop calendar, fertilizer, and soil, which were harmonized at the same resolution (see table S1 available online at stacks.iop.org/ERL/17/024028/mmedia). The grid level simulated yields were aggregated to country level and compared with country level observed yields (text S1). The standardized root mean square error (RSR) criterion proposed by Singh et al [49] was used as the objective function to evaluate the model performance when comparing country-level simulated yield with country-level yield from the Food and Agriculture Organization (FAO) [50] over a 33 years period (1980–2012). Further information on EPIC crop-related processes is available in Williams et al [45]; the procedure for calibration of EPIC⁺ is found in supplementary text S1 and in Kamali et al [20].

Different drought indices offer different information about a given drought phenomenon. A meteorological drought resulting from a precipitation deficit may develop rapidly, while a deficit in soil moisture as a measure of agricultural drought may occur with a time lag. Recently, Rigden et al [51] found that root-zone soil moisture is a better proxy than using standard approaches based on temperature or precipitation for agricultural purposes and yield prediction. In this study, we used and compared two standardized indices—the standardized precipitation index ((SPI) [52]) and the standardized soil moisture index (SSI [53, 54])—to justify the suitability of soil moisture for our analysis in the study area and also to quantify both variants and describe their severity.

The indices can be calculated at different time scales. As also suggested by Kamali et al [20], the 3 and 6 months time scale were identified to be most suitable for SSA. We here also selected these two time scales depending on the highest correlation with crop yield. Table 1 shows the six levels of drought severity grouped according to the World Meteorological Organization standards [55] used for characterizing SPI and SSI.

Table 1. Six categories of the SPI, the SSI, the CSI (crop sensitivity index) and YRisk (risk of yield loss under drought).

Ranges	SPI, SSI	CSI	YRisk
1.5–3.0	Extremely wet	Very high yield	Very low
1.00–1.49	Very wet	High yield	Low
0–0.99	Moderately wet	Moderately high yield	Moderately low
$-$0.99–0	Moderate drought	Moderately low yield	Moderately high
$-$1.49 to $-$1.00	Severe drought	Severely low yield	High
$-$3.0 to $-$1.5	Extreme drought	Extremely low yield	Very high

The CSI for each year t was calculated using the yield values (${Y_{t,i}}$) simulated in each grid cell (i) for each crop as:

$$\begin{equation}{\text{CS}}{{\text{I}}_{t,i}}\ =\ \frac{{{Y_{t,i}} - \bar Y}}{{{\sigma _Y}}},\end{equation} \tag{ 1 }$$

where $\bar Y\,{\text{ and }}{\sigma _Y}{ }$ indicate the mean and the standard deviation of the yield values, respectively over 33 years (1980–2012) for each crop and each grid point, the ${\text{CSI}}$ time series has a mean of 0 and a standard deviation of 1. Generally, EPIC estimates different stress factors caused by water, nutrients, temperature, aeration, and radiation [45]. These stress factors vary between 0 and 1 (with 1 indicating no stress and 0 the highest stress). The potential biomass and yield will be constrained by using the lowest factor among the stresses. The impacts of nitrogen stress were of course considered during model calibration. In order to explore the impact of drought on crop yield (crop drought risk analysis), we deployed the calibrated crop model, and simulated three crops across grids in SSA under a condition in which the model's response to nitrogen stress was turned off. Simulating yields under no nitrogen stress provided the possibility to directly assess the impact of drought and high-temperature stresses in different regions. The simulated yields affected by water stress only were then used for the probabilistic modeling of crop drought risk. Similar to SPI and SSI, six categories of yield levels were defined to characterize and classify CSI (table 1).

We modeled the risk of yield loss under drought (YRisk) probabilistically using Copula functions. Copula functions are multivariate cumulative distribution functions used to describe the dependence between random variables. The aim of this step is to estimate the responses of crop yields to different droughts of various severity levels. The main advantage of Copula-based methods is that they enable researchers to model the dependence structure between variables with different, including non-normal, marginal distributions [56]. We therefore described the relationship between a drought index (here, the SPI or the SSI) and the CSI using a Copula function [56] (figure 1). Assuming that X is the random variable of drought derived from soil moisture deficit or precipitation deficit, and Y is the random variable of crop yield levels, the joint probability distribution between the two (X and Y) can be expressed as:

$$\begin{equation}{\text{Prob}}\left( {X \le x,{ }Y \le y} \right) = {f_{X,Y}},\end{equation} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

where ${f_{X,Y}}$ denotes the joint probability distribution of the soil moisture (or precipitation) level and the crop yield level, we used a nonparametric joint distribution concept to construct the joint probability distribution. In doing so, the empirical joint probability in the bivariate case using the Gringorten plotting position formula [57] emerged as:

$$\begin{equation}{f_{X,Y}} = \frac{{{m_k} - 0.44}}{{n + 0.12}}\,\end{equation} \tag{ 3 }$$

where n is the number of observations (here, 33 years) and ${m_k}$ is the number of occurrences of the pair $\left( {{x_k},{y_k}} \right)$ for ${x_i} \le {x_k}$ and ${y_i} \le {y_k}$ [58]. Generally, variables such as temperatures, precipitation, and soil moisture can be non-stationary especially under a changing climate [57]. This makes it difficult to interpret fitted extreme value of distribution functions using parametric methods. Therefore, the application of nonparametric methods is appropriate [54]. The main advantage of the nonparametric methods is that they do not require making presumptions about the distribution of the variables. Besides, they alleviate the computational burdens in fitting the parametric distribution. Having determined the joint distribution, the risk of crop yield loss given a drought event (YRisk) is defined as [54]:

$$\begin{equation}{\text{YRisk}} = {\emptyset ^{ - 1}}\left( {{f_{X,Y}}} \right),\end{equation} \tag{ 4 }$$

where $\emptyset$ is the standard normal distribution function. Almost all values of ${\text{YRisk}}$ vary between −3 and 3 where −3 indicates the highest risk and 3 indicates the lowest risk. YRisk was then grouped into six categories of severity showing the risk level of crop loss under a drought event (table 1).

We quantified how much each of the various drivers (climatic conditions and soil types) contributed to the risk of crop-yield loss (YRisk) by computing variance-based sensitivity indices [14, 59]. We used variance-based methods in which the variance of output (here YRisk) is decomposed into fractions and are attributed to inputs or sets of inputs (to different factors). This allows us to compare the levels of influencing factors on each crop. The variance based method can deal with nonlinear responses, and they can measure the effect of interactions in non-additive system. The decomposition analysis is conducted considering soil type, precipitation amount, and temperature during the main development stages of the crops. Accordingly, the contributing factors included: (a) the soil type; (b) precipitation from sowing to emergence (PCP_SE); (c) precipitation from emergence to flowering (PCP_EF); (d) precipitation from flowering to grain filling (PCP_FG); (e) the average temperature from sowing to emergence (TMP_SE); (f) the average temperature from emergence to flowering (TMP_EF); and (g) the average temperature from flowering to grain filling (TMP_FG). Specifically, two sensitivity indices,—the first-order or main effect (E_i ) and the total effect (TE_i )—were applied to estimate the variance caused by one factor from the variance caused by the interaction of the rest of the factors:

$$\begin{equation}{E_i} = \frac{{\operatorname{var} (E[{\text{YRisk}}|{V_i}])}}{{\operatorname{var} \left( {{\text{YRisk}}} \right)}}\end{equation} \tag{ 5 }$$

and,

$$\begin{equation}{\text{T}}{{\text{E}}_i} = 1 - { }\frac{{\operatorname{var} \left( {E[{\text{YRisk}}|{V_{i - 1}}]} \right)}}{{\operatorname{var} \left( {{\text{YRisk}}} \right)}},\end{equation} \tag{ 6 }$$

where ${\text{var}}$ is the variance and $E[{\text{YRisk}}|{V_i}]$ denotes the expected value of YRisk across factor ${V_i}$, while $E[{\text{YRisk}}|{V_{i - 1}}]$ is the expected value of YRisk across all factors except ${V_i}$. The main effect E_i measures the effect of varying factor alone. On the other hand, not only does the TE_i include the variance of an input, it also accounts for the variance created by interaction with other parameters. The total effect index is a summary sensitivity measure inclusive of interactions effect of any order [60].

The yields of maize, millet, and sorghum were simulated by EPIC⁺ at grid level. The simulated yields were then aggregated at the country scale and compared with the country-level yields reported by the FAO [50] over the 1980–2012 period. The goodness of fit (R²) values of 0.94, 0.90, and 0.96 for maize, sorghum, and millet, respectively with their corresponding p-value close to 0.0 for all crops show the level of agreement between the country-level simulations and observations (figures 2(a)–(c)). However, aggregated results can mask the regional uncertainties of the model simulations as well as their temporal trends. We then analyzed the temporal dynamics of simulated yields during 1980–2012 and their uncertainty using the SUFI-2 algorithm at the country level, relative to observations. For some regions (South Africa, Kenya, Benin, and Chad) as shown in figure S1, the RSR values are smaller than 2, indicating acceptable yield simulations. The p-factors with values around 0.5 or more in all countries indicate that over 50% of the data is bracketed within the 95PPU band, which is acceptable given the complex nature of the problem (table 2). Furthermore, the r-factor values, representing the width of the 95PPU band, are smaller than 2.5 (table 2), indicating reasonable levels of uncertainty for crop simulations [20].

Zoom In Zoom Out Reset image size

Using the SPI and the SSI, we identified the drought events that affected most parts of SSA during 1980–2012. Both indices were consistent with each other, revealing 1983, 1984, 1987, 1992, and 1995 as the drought years. However, there is a slight discrepancy in terms of severity levels (supplementary figure S2). As shown in figure S2, the SSI showed with a higher severity than the SPI for all five drought events. The extreme drought event of 1983 affected all of SSA, but with less severity in eastern Africa. A year later (1984), eastern African countries experienced extreme drought, while the drought in western Africa was less severe than in 1983. The drought of 1987 influenced mostly southern and western Africa. The 1992 drought developed with higher severity in the SSI as compared with the SPI, affecting most regions of SSA particularly with extreme severity in the southern and eastern parts. Another drought event was identified in 1995, mostly affecting central SSA (figure S2).

Next, we looked at the relationship between the drought indices (SSI or SPI) and CSI at the country level for the five identified drought years (figures 3(Aa–o)) and five other randomly selected years of 1996, 1999, 2000, 2006, and 2012 (figures 3(Ba–o)). As shown in figures 3(Aa–o) for five drought years, the relationships between the SSI and the CSI indicate a consistent decrease in crop yields as a result of decreases in SSI values. Overall, the CSI variability could be explained by the SSI drought index to a large extent in all five selected drought years. However, depending on the crop and country, the CSI responded differently to various severities of the drought events. Generally, the correlation coefficient between two indices (r) exhibited higher values for maize and millet than for sorghum, for which the (r) values were lower than 0.5 in all five drought events. The high correlation means that a drought of a specific severity results in a corresponding yield loss of the same severity. However, lower r values indicate that yield losses may also be related to factors other than drought.

Zoom In Zoom Out Reset image size

In 1983, for example, the correlation coefficient values were 0.71 for maize and 0.70 for millet, whereas it was only 0.41 for sorghum. Analyzing the relationship between the SSI and the CSI at the country level also indicated different responses from one region to the other. For instance, in 1984 eastern African countries exhibited higher CSI values under the influence of the SSI (figures 3(Ag) and (Ai)), meaning that the drought that year resulted in less yield reduction. A similar pattern emerged for the western African countries in the same year. This implies that the occurrence of a drought event of similar severity may not necessarily lead to the same magnitude of yield reduction. The relationship between the SPI and the CSI during these drought years shows comparable trends, but with different r values (supplementary figure S3). We further, investigated the feasibility of the proposed approach in five randomly selected years (i.e. 1996, 1999, 2000, 2006, and 2012) throughout the study period (figures 3(Ba–o) and supplementary figures S3(Ba–o)). Comparing drought years with non-drought years show that both SPI and SSI can clearly distinguish the events from each other. It is also the case for the regional drought events, for instance, in the year 2000 some eastern African countries were experiencing severe to extreme droughts with high to very high YRisk while other African countries were in a better condition in terms of SSI and SPI with low to very low YRisk level (figures 3 and S3 ((Bc), (Bh), (Bm))). However, under moderate drought to moderately wet conditions (SPI or SSI between −0.99 and 0.99), CSI varies in a wider range and the correlation coefficient (r) between drought indices and CSI are low, again implying the importance of other factors in the yield production during these conditions.

We then looked at YRisk (equation 4) at the grid cell level for the whole SSA. The results for the drought years (figure 4(A)) and randomly selected years (figure 4(B)) showed that YRisk varied substantially by location, crop, and year. For example, the drought of 1984 resulted in a high risk of sorghum loss in almost all western African countries, while the risk of maize yield loss was higher for the Sahel countries (of western Africa) than for the other western African countries. Similar differences among different crop drought risks were seen in 1983, where YRisk for millet was slightly lower. By visually comparing the spatial distribution of risk severity (figure 4(A)) with drought severity (figure S2), a strong correspondence between the occurrences of drought and yield reduction emerged. This is also found in the case of regional droughts, for instance in the year 2006, the western countries were experienced drought with high to very high YRisk (figures 3 and 4((Bd), (Bi), (Bn))). However, the levels of severity varied remarkably in both cases, confirming the appropriateness of applying spatially distributed and probabilistic approaches in analyzing the drought effects on crops.

Zoom In Zoom Out Reset image size

The shape of the kernel density functions was compared for the three crops (figure 5), exhibiting variable shapes depending on the crop and country. As shown in figure 5, for the sorghum crop there is a shift in the density function to the left side, revealing a higher risk of yield loss due to drought compared to other crops. The higher risk of drought for sorghum (compared to the other two crops) was more evident in countries like the Central African Republic (Central Africa), Gambia, Ghana, and Somalia than in the rest of SSA. In most countries, we also see elongated tails on the left side compared to the right side, which confirms the higher probability of YRisk with higher density.

Zoom In Zoom Out Reset image size

Next, we conducted a decomposition analysis based on the main and total effect indices for the three crops during drought years to understand which of the seven factors considered here contributed the most to the risk of yield loss under drought conditions. The results show that factors contributing to the level of risk were the same for all three crops, but their relative contributions varied slightly from one crop to another (figure 6). Generally, the amount of precipitation during vegetation and grain filling, i.e. PCP_EF and PCP_FG, contributed the most to YRisk. The relative contribution of temperature during vegetation and grain filling, i.e. TMP_EF and TMP_FG, were highest for sorghum, followed by millet and maize, respectively. The results for the total effect showed higher interactions of soil type with precipitation and temperature variables (figure 6).

Zoom In Zoom Out Reset image size