Drivers of Subseasonal Forecast Errors of the East African Short Rains

The ‘short rains’ in East Africa from October to December have significant year‐to‐year variability. Their abundance or deficiency is often associated with floods or droughts for which early warning is crucial, though even in normal seasons skillful forecasts facilitate planning and preparedness. Here we study the relationship between initial‐state sea surface temperatures and subseasonal rainfall forecast errors in the European Center for Medium‐Range Weather Forecasts model in the region. We demonstrate that the initial mode of the Indian Ocean Dipole (IOD) is a partial control on the rainfall error in weeks 3–4. This relationship is also clear on the seasonal scale, exemplified by too‐wet forecasts during the 2015 season when the IOD was positive, and too‐dry forecasts in 2010 when it was negative. Our results provide an entry point for model improvement, and we show that a priori forecast corrections based on the initial IOD index are feasible.

routinely disseminated in East Africa. This means that following the pre-season outlook, the next rainfall forecast update is often issued only a few days ahead of the hazard (e.g., Kilavi et al., 2018), limiting the potential for risk management.
The lack of subseasonal forecast provision is unfortunate, as compared to most of the globe East Africa is a hotspot of forecast skill (de Andrade et al., 2020;Vigaud et al., 2019). The strong predictability arises in large part from the Madden-Julian Oscillation (MJO), which shows a strong teleconnection with regional rainfall: models with skillful MJO predictions also capture the teleconnection to rainfall over East Africa (MacLeod et al., 2021).
An impediment to the uptake of subseasonal forecasts is a lack of robust understanding of model performance in the region. Recent studies are addressing parts of this knowledge gap by demonstrating levels of model skill and evaluating model representation of teleconnections (Hirons & Turner, 2018;King et al., 2021;Walker et al., 2019). However, most studies focus on unconditional model errors (e.g., mean model biases) rather than conditional model errors, that is, forecast errors associated with particular climate states (e.g., Kolstad et al., 2020;Matsueda & Palmer, 2018;Minami & Takaya, 2020). Understanding of conditional model errors can give two benefits supporting the reliability of operational forecasts. First, it can make possible a context-specific a priori error correction: knowing that the model will behave a certain way under a certain climate state allows a more intelligent use of model output. Second, it provides targeted error feedback to model developers, which can potentially lead to improvement of simulation realism in the region and higher model skill.
Here, we study the conditional model error of subseasonal forecasts of the East African short rains. In particular, we investigate the relationship between the initial state of the Indian Ocean Dipole (IOD; Marchant et al., 2007) and subsequent forecast error developments. Our aim is to enhance understanding of model performance, which we hope may contribute to raising the level of trust in forecasts on decision-relevant time scales.

Data
The European Center for Medium-Range Weather Forecasts (ECMWF) monthly forecasts are issued twice a week. The model version used here is CY46R1, which has a native horizontal atmospheric grid spacing after day 15 of about 32 km, and the ocean model is NEMO3.4.1, with a native 0.25° grid spacing. Although the 32 km grid spacing is insufficient for resolving small-scale convective systems, it is sufficient for reproducing the large-scale rainfall characteristics driven by ENSO and the Walker Circulation over the Indian Ocean. A set of 11-member hindcasts is produced for each forecast, initialized for the same day and month of the previous 20 years. From the S2S database (Vitart et al., 2017), we use the ensemble means of the first 11 members of the (51-member) forecasts for which all lead times between 15 and 28 days (3-4 weeks) occurred between October 1 and December 31, 2019, as well as their associated hindcasts. This comprises 23 × 21 = 483 model runs. The variables studied are precipitation (hereafter referred to as "rainfall"), zonal wind at 850 and 200 hPa, and SST, all interpolated to 1.5° grids, which is how the atmospheric data are stored in the S2S database. For the same two-week periods, we also calculate accumulated rainfall from the Climate hazards infrared precipitation with stations (CHIRPS) data set (Funk et al., 2015), interpolated to the 1.5° model grid. For the other variables, we use ERA5 reanalysis data (Hersbach et al., 2018a;2018b; as our reference.

Standardization
To account for spatial heterogeneity and model bias, we post-process the rainfall data in several steps. First, we standardize the 483 rainfall time series for each grid point from both the forecasts and the CHIRPS data by subtracting the datewise climatological mean and dividing by the datewise interannual standard deviation. The smooth the climatologies, they are based on data from all (171, on average) dates up to 14 ordinal days before or after the ordinal day of the reference date. Second, we calculate the area-average of the standardized forecast anomalies within our focus region -between 30°E and 52°E and between 10°S KOLSTAD ET AL.

SST Indices
SST data from ERA5 are used to calculate an IOD index (I henceforth) in two steps. First, we compute area-averaged standardized SST anomalies (using the same climatological smoothing as for the rainfall data) in the western (50°E to 70°E and 10°S to 10°N) and eastern (90°E to 110°E and 10°S to 0°) parts of the Indian Ocean. Second, we take the difference between the western and eastern time series, and third, we standardize the resulting time series. Although our main focus is on the IOD index, we also calculate a NINO3.4 index as area-averaged standardized SST anomalies between 170°W and 120°W and between 5°S and 5°N. The boundaries of the SST regions are shown in Figure 2a.

Significance Testing
To account for intraseasonal autocorrelation, we use bootstrapping to create sets of 1,000 artificial time series to contrast with the actual time series. Each time series consists of data from 23 initial dates per year, and in the artificial time series, we replace each set of 23 elements with a set from a random year between 1999 and 2019 (with replacement). This retains the intraseasonal autocorrelation. The null hypothesis that actual correlation or regression coefficients could just as well have resulted from the randomized, artificial sets can be rejected at the 5% significance level (which is used throughout) if the actual coefficient falls outside the interval between the 2.5th and 97.5th percentiles of the artificial set.

Linear Prediction Model
In Section 3.2, we use a linear regression model to predict O based on the predictors F and I: (1) The coefficients are estimated (by minimizing the residual  using ordinary least squares) for each forecast time by using F and I for all the other forecast times ("leave-one-out" cross validation), and we call the best estimates ˆF a and ˆI a , so that the prediction of O can be written as: Note that we do not compute out-of-sample climatologies for each element of the time series, as the high number of data points these are based on (171, on average) makes the in-sample effect negligible.

The Short Rains
We first present an overview of the short rains' relative importance in the annual rainfall cycle, as well as its representation by the model. The ratios in Figure 1a were obtained by dividing the average rainfall during all the 483 OND forecast periods by the average rainfall during all similar periods throughout the year. The OND rainy season is most dominant compared to the rest of the year in Somalia, southern Ethiopia, and eastern Kenya, but ratios above 1 also occur in other parts of the study region, including its southwestern part. 10.1029/2021GL093292 3 of 10 Figure 1b shows that the anomaly correlation between the forecasts and the regridded CHIRPS rainfall is significant in most of the whole study region. The minimum, maximum and median correlation coefficients within the region are -0.01, 0.61, and 0.31, respectively. On the spatially aggregated scale, the (significant) correlation between F and O is 0.55 (Figure 1c). This is substantially higher than the median of the grid point correlations, and it means that the ensemble mean forecast explains 30% of the rainfall variability 3-4 weeks in advance on the aggregated scale. The best fit line is less steep than the identity line (i.e., where  x y), indicating that F needs calibration. Figure 1d shows a scatter plot of I versus O. The correlation is 0.51, which is lower than the correlation between F and O. In 2019, the high I values are associated with a wide range of O values from about zero to the five largest values in the study period. There is clearly information in I not captured by F, suggesting that a prediction model based on combining F and I might be worthwhile to test (see Section 3.2).
Although not shown, we note here that the correlation between O and the initial NINO3.4 index is only 0.23. As a result, our focus hereafter is on the IOD index.

Linking Initial SSTs and Rainfall Forecast Errors
While it is well-known that East African OND rainfall is correlated with both ENSO and the IOD, Figure 2a demonstrates that the forecast error  is also lag-correlated with initial tropical SSTs. This means that the ECMWF model has a systematic bias conditional on the initial state, and it also suggests that initial SSTs can potentially be used to correct the model forecasts a priori.
To test whether the forecasts can be corrected based on I, we assess the performance of Ô , the linear prediction of O based on F and I (see Section 2.5). A scatter plot O versus Ô is shown in Figure 2b. The correlation between O and Ô is 0.58, which is significantly higher than the correlation between O and F (0.55; see Figure 1c) according to a bootstrapping test (Section 2.4). This means that including I as a predictor in combination with F is indeed worthwhile.
As the correlation between I and  is relatively high (0.24) and significant (see also Figure 2a), we were curious about why Ô does not represent an even larger improvement with respect to F. It could be related to nonlinear aspects of the relationship between  and I. Figure 2c shows boxplots (e.g., Krzywinski & Altman, 2014) of  for four bins organized by ascending I values (the leftmost bin is based on the forecasts initialized when I was in its lower quartile, and so on). The linear part of the relationship between I and  is reflected by the increasing median values from left to right, but it is also clear that the variance of  within  We also tried including the initial NINO3.4 index as a predictor of O, alongside F and I, but its regression coefficient was nonsignificant. This is commensurate with the low correlation between O and the NINO3.4 index (Section 3.1).

Drivers of Rainfall Forecast Errors
To better understand the dynamics behind the rainfall forecast errors, we show in Figure 3a Figure 3c shows negative correlations between zonal wind errors at 200 hPa over the continent and off the East African coast, and positive correlations over the eastern Indian ocean. This is an indicator of a surplus of divergence aloft in the western part of the Indian Ocean, which is linked to too strong (weak) convection below 200 hPa in the model when  is positive (negative).
We now investigate the relationship between I and the SST and zonal wind forecast errors by comparing the correlations shown in the bottom row of Figure 3 with the correlations in the top row. We emphasize that  does not enter into the calculations of the correlations in the bottom row. Still, the correlation between the SST errors and I (Figure 3d) in the northern Indian Ocean agrees well with the correlations between the same SST errors and  (Figure 3a). For 850-hPa zonal wind, shown in Figure 3e, there is also good KOLSTAD ET AL.   (Figure 3f) is strong over the Indian Ocean, the magnitude of the negative correlation with  over land is not reproduced for I aloft. This indicates that model errors in representing the variability of the convergence and divergence over the African continent are a determinant of subseasonal forecast errors over East Africa, and that the origin of these errors is only weakly related to the initial state of the IOD.

Seasonal Means
The aforementioned intraseasonal autocorrelation of forecast errors (e.g., Figure 1c), in combination with a significant link between the initial IOD index and  , suggest a possibility for a priori correction of subseasonal forecast error at the start of the season (or even before, given skillful long-lead SST forecasts). Although based on only 21 seasons, the correlation between the seasonal means of  and I, shown in Figure 4a

Summary and Discussion
We have demonstrated a robust link between subseasonal rainfall forecast errors over East Africa, tropical SSTs, and the atmospheric circulation. A symmetry has been identified, where positive (negative) IOD states in the initial conditions are associated with too-strong positive (negative) rainfall anomalies over East Africa in the model after 3-4 weeks.
KOLSTAD ET AL.
10.1029/2021GL093292 6 of 10 A linear correlation approach was used to link the rainfall forecast errors to forecast errors of oceanic and atmospheric drivers (see Figure 3, top row). Focusing on the too-wet forecasts (the too-dry forecasts mainly show "opposite" behavior), these were found to be associated with concurrent model errors in the western Indian Ocean region, specifically in the form of positive SST errors, coupled with too strong divergence aloft and too strong easterly zonal wind anomalies at low levels. Furthermore, we found that the too-wet forecasts are associated with positive low-level zonal wind errors stretching into East Africa from the west. Previous studies have repeatedly demonstrated an association between easterlies over the Indian Ocean and enhanced rainfall and flooding in East Africa (Black et al., 2003;Hastenrath, 2007;Nicholson, 2017). However, the results presented here are unique by linking this state not to rainfall totals, but to errors in the rainfall forecast.
Another key finding is that the correlations between the errors in the forecasts of SSTs and zonal winds and the rainfall forecast error exhibit many similarities to the correlations between the initial IOD index and same variables' forecast errors (Figure 3, bottom row). This suggests oceanic and atmospheric pathways linking the initial IOD index with subsequent rainfall forecast errors. The forecast errors during the two OND seasons with the highest (in 2019) and lowest (in 2010) mean IOD index in the study period confirmed a linkage between the IOD index, rainfall forecast errors and forecast errors of SSTs and zonal winds on the seasonal scale.
KOLSTAD ET AL.

Too-dry forecasts Too-wet forecasts
These results have two important consequences. The first is that they suggest an opportunity to make a state-dependent correction to the forecasts. This could be done in a systematic way whereby forecasts initialized in anomalous SST states are calibrated, with a correction to the ensemble mean rainfall forecast. As a proof-of-concept exercise, this was done here with a linear adjustment of the forecast based on the initial IOD index in Section 3.2, and the corrected forecast had significantly higher skill than the uncorrected forecast. We emphasize that our intention is to demonstrate the potential for using initial SSTs to create a hybrid empirical/dynamical forecast model, and not to optimize such a model. Thus, we did not attempt to correct the rainfall forecasts on the local scale (see skill in Figure 1b), but the promising results for the aggregated geographical scale suggest that it might be possible. But even without a systematic calibration, knowledge of the systematic relationship between initial SSTs and rainfall forecast errors is likely to be useful to experienced forecasters, who may incorporate this into their subjective interpretation when developing forecast products.
The second main consequence is that our results provide modelers with a point of entry for addressing systematic model errors, which would offer hope of improved simulation realism and downstream delivery of more societally relevant model forecasts in the region. An investigation might proceed by evaluating the presence of the error in previous model versions to see if any recent developments have improved (or exacerbated) the issue. Contrasting with the performance of other models may also offer insight. Following this a smaller test case could be identified, such as the forecasts for the particular seasons discussed here. Reforecast experiments could be run to test the impact of changes in model setup, such as resolution, convective parameterization or atmosphere-ocean coupling parameters. Alongside this a close inspection of the day-to-day evolution of the ocean-atmosphere state during particularly extreme errors may provide insight into their origin.

Data Availability Statement
The CHIRPS data (Funk et al., 2015) are available from the Climate Hazards Center at the University of Santa Barbara, USA at https://www.chc.ucsb.edu/data/chirps. The ERA5 data are available from the Copernicus Climate Data Store (Hersbach et al., 2018a;2018b). This work is based on S2S data. S2S is a joint initiative of the World Weather Research Programme (WWRP) and the World Climate Research Programme (WCRP). The original S2S database (Vitart et al., 2017) is hosted at ECMWF as an extension of the TIGGE database at https://apps.ecmwf.int/datasets/