Machine Learning Prediction of the Rainy-Season Onset in Southeast Africa

1. Introduction

The onset of the rainy season is a key climate event in Southeast Africa with implications for rainfed agriculture, water management, power production and disaster preparedness. Forecast information is particularly valuable when it indicates the likely timing of seasonal rains which is crucial for mitigating the adverse impacts of delayed or erratic rainfall. However, predicting the onset of the rainy season remains a challenging task due to the complex interplay of atmospheric, oceanic, and land-surface processes that influence regional climate patterns (Roy et al., 2023; Palmer et al., 2023; Roy and Troccoli, 2024).

In recent years, advancements in machine learning and data-driven approaches have opened new possibilities for improving the skill of seasonal forecasts. The Random Forest (RF) method, a widely used machine learning algorithm (Breiman, 2001), has shown great potential in handling the complexity of climate data, especially when trained on large amount of data, as demonstrated by Gibson et al. (2021), where RF achieved skill on par with dynamical models and offered interpretable insights into physical drivers.

This study develops gridded forecasts of the October-November-December (OND) rainy season in Southeast Africa, covering the countries of Tanzania, Malawi, Mozambique, Rwanda, and Zimbabwe, with a focus on Tanzania, through close collaboration with the Tanzania Meteorological Authority (TMA), as part of the 4-year EU-funded FOCUS-Africa project. These countries were selected because they formed the core geographical focus of the FOCUS-Africa activities, and a regular rectangular study domain was adopted to encompass them consistently. While the analysis presented here focuses on this domain, the methodology is generic and can readily be extended to neighbouring countries and regions, including Kenya and Ethiopia. TMA’s local expertise significantly informed both current methods for onset calculation and stakeholder requirements for forecast applications.

Two internally consistent data systems were constructed: one based on Multi-Source Weighted-Ensemble Precipitation (MSWEP) (Beck et al., 2019) and one based on Climate Hazards Group InfraRed Precipitation with Stations (CHIRPS) (Funk et al., 2015). For consistency, the same precipitation dataset was used to derive both the onset target and the associated lagged-rainfall predictors. These local inputs were combined with observed and forecast ENSO and IOD indices and with the dynamical seasonal onset forecasts. ENSO and the IOD were included because they are major large-scale drivers of OND rainfall variability and onset timing in Southeast Africa. ENSO can alter the Walker circulation and regional moisture transport, while the east–west sea-surface-temperature gradient associated with the IOD influences convection over the western Indian Ocean and East Africa. These relationships are season-, location-, and lead-time-dependent rather than deterministic (Wenhaji Ndomeni et al., 2018; Kolstad and MacLeod, 2022; Roy and Troccoli, 2024). Observed indices describe the large-scale state at forecast issue time, whereas forecast indices describe its expected evolution during the lead-up to OND. Their inclusion is therefore physically motivated but does not imply a causal interpretation of RF predictor importance.

Building on previous work (Zampieri et al., 2023; Roy et al., 2023), this study evaluates whether statistical post-processing can add useful information to dynamical seasonal onset forecasts across different initialization months.

The accumulated-anomalies method, proposed by Zampieri et al. 2023, is used here to create a more generic onset metric that can be applied across larger regions and different domains than the threshold-based method currently used by local institutions such as TMA. By comparing the model’s predictions for the 2023 and 2024 rainy seasons against ECMWF system 5.1 (s5.1) dynamical forecast systems, the study demonstrates the potential of the RF method to support decision-making in climate-sensitive sectors across Southeast Africa. ECMWF s5.1 was chosen as the primary baseline because multiple independent evaluations have shown that, among Copernicus seasonal forecast systems, ECMWF generally provides the highest precipitation forecast skill for East and Southern Africa, outperforming DWD and CMCC in terms of monthly anomaly correlations, predictive capability for the OND rains, and operational use across the region (Gebrechorkos et al., 2022; Tefera et al., 2025; Endris et al., 2021). As our RF model incorporates dynamical-model onset forecasts, including ECMWF-derived information, this comparison should be interpreted as an assessment of statistical post-processing relative to a raw dynamical forecast, rather than as a comparison between two fully independent forecasting systems.

2. Data and Methods

Two precipitation products were used to construct parallel historical datasets of the OND rainy-season onset: MSWEP and CHIRPS. The two products were not merged into a single training dataset. Instead, each product defined an internally consistent system: onset dates were derived from that product and its corresponding antecedent-rainfall predictors were calculated from the same product. Together, these provided two independent training datasets based on robust daily precipitation estimates. A comprehensive assessment of available precipitation datasets was carried out as part of FOCUS-Africa, identifying these as preferred options for Tanzania, owing to their long, spatially complete precipitation records and suitability for gridded onset estimation (Mliwa et al., 2025, in press). The characteristics and preprocessing of each dataset are described below.

MSWEP provides global precipitation estimates with a 0.1° resolution and three-hourly temporal resolution, combining satellite data, ground station measurements, and reanalysis products. It offers high temporal frequency (3-hourly data), making it ideal for hydrological applications (Beck et al., 2019). It was used both to provide MSWEP-based onset dates and to calculate the corresponding lagged-rainfall predictors.

CHIRPS, a high-resolution (0.05° grid) dataset, is designed for drought monitoring and combines satellite-derived estimates with ground station data. It provides daily rainfall estimates, particularly useful for regions with sparse observational networks like Southeast Africa (Funk et al., 2015). For this study, the regridded version of CHIRPS available at 0.25° resolution was chosen for computational purposes. The MSWEP data were then bilinear interpolated onto the CHIRPS grid, resulting in a common 0.25° grid spanning 27.875° S–0.125° S and 24.125° E–40.875° E. Dynamical seasonal forecast inputs were subsequently interpolated onto the same grid, with the same method, before model training and evaluation.

For each precipitation product, historical OND rainy season onset dates were calculated independently for every grid point and year. These were calculated as proposed by Zampieri et al. (2023) by identifying the point at which accumulated daily precipitation anomalies, calculated from the first of July, reached their minimum before a sustained period of above-average rainfall began. The climatology was defined as the mean daily rainfall from 1 July to 15 January, with reference to the 1993-2016 hindcast period. Formally the daily anomaly was calculated as:

α_t  =  P_t  −  P_c

(1)

Where P_t is daily precipitation and P_c is the climatological daily mean. The accumulated anomaly from 1st of July was then calculated as:

$$A_t = {\sum }_{s=1}^t{a}_s$$

(2)

Where A_t is the accumulated precipitation anomaly on day t and as is the daily precipitation anomaly on each day (s) from 1st of July to day t.

To account for the bimodal rainfall regime present in parts of the study area, particularly in northern Tanzania, the method was adapted. In these zones, an initial increase in accumulated daily rainfall anomalies was accepted as the onset if it resulted in a local maximum that exceeded the subsequent minimum by more than 40 mm, thus recognising the occurrence of an early, distinct rainy season, even if another rainfall minimum followed. The 40 mm separation criterion was part of the locally informed onset definition adopted for this study and was applied consistently at every grid point and for every year, rather than fitted separately by location or forecast start month. This approach identifies a distinct early wet period rather than automatically assigning onset to the later minimum. This method was chosen instead of onset estimation methods based on fixed daily-precipitation thresholds, as such methods are typically tailored to specific countries or areas and might not perform consistently across domains as large as the one considered here. Figure 1 shows an example of the onset calculation for the year 1997 from a randomly selected grid point.

Our RF model used 23 predictors, drawn from two main sources: observational data and dynamical forecast data (Table 1). With the data being retrieved from Copernicus Data Store (CDS), National Oceanic and Atmospheric Administration (NOAA), International Research Institute for Climate and Society (IRI) and GloH2O.

The observational predictors were restricted to information available before each forecast initialization date. They include ENSO 3.4 index (sea surface temperature, SST, anomalies in the central Pacific) and the Indian Ocean Dipole (IOD) index, which reflects SST anomalies across the Indian Ocean. The ENSO 3.4 was used with a lags of 1-3 months prior to forecast initialization and IOD with lags of 1-4 months. Forecast ENSO and IOD predictors represent information available at the forecast issue date about the expected evolution of large-scale sea-surface-temperature conditions before and during the OND season. For a forecast initialised in month m, fc_L1, fc_L2, and fc_L3 denote the ensemble-mean index forecasts for the first, second, and third months following that initialisation, respectively. Observed and forecast indices have complementary roles: observed lags describe the realised large-scale climate state before forecast initialisation, whereas forecast indices describe the modelled evolution of that state into the forthcoming months. Since they may be correlated, their inclusion does not assume that they provide independent information. Rather, our RF is allowed to use either or both sources of information where this improves out-of-sample prediction; their importance is not interpreted as a causal effect. In addition, local precipitation predictors derived from MSWEP or CHIRPS with lags of 1-3 months were included to capture local rainfall variability. The rainfall predictors followed the parallel-data-system approach: MSWEP monthly averages were used when the onset was calculated from MSWEP daily accumulated anomalies and CHIRPS was used similarly for CHIRPS-based onset estimates.

The forecast data used as predictors in our RF model consisted of onset forecasts from three distinct dynamical seasonal forecast systems: the European Centre for Medium-Range Weather Forecasts (ECMWF) system 5.1 (s5.1), the Deutscher Wetterdienst (DWD) s2.1, and the Centro Euro-Mediterraneo sui Cambiamenti Climatici (CMCC) s3.5. For each forecast system, the predicted onset of the rainy season was derived using the same accumulated rainfall anomaly method applied to the observational target. Forecasts from individual ensemble members were aggregated by calculating the mean predicted onset date, which was then used as a predictor in our RF model. Figure 2 shows an example forecast of accumulated daily anomalies from the ECMWF s5.1, with predicted onset dates from individual members marked by green circles and the ensemble-mean onset, used as a predictor, indicated as a dashed vertical line. Combined, these predictors allowed our RF to function as a statistical post-processing system that augments dynamical model forecasts using observational information, rather than as an entirely independent forecast system.

In addition to the forecasts of onset, forecasts of ENSO3.4 and IOD were also generated from SST forecasts from the required domains. These indices were calculated for lead times of 1-5 months for all initialization dates, first by averaging across all ensemble members for each model, and then by averaging across models.

The RF model was chosen as the ML model for this study, after comparison with Bayesian Ridge and multilayer perceptron models under the common-predictor leave-one-year-out hindcast framework; it was not assumed a priori to be the most skilful method in every setting. Furthermore, it was chosen due to its ability to handle complex datasets with numerous predictors. The RF model generates multiple decision trees during training, each providing a prediction. The final output is an ensemble of all individual tree predictions, ensuring a probabilistic forecast that reflects its uncertainty and can capture non-linear relationships between the predictors.

Our RF model was trained at each of the 7616 grid points over the whole study domain. Separate RF models were trained at each grid point. Each model contained 250 trees and used squared-error reduction as the tree-splitting criterion. Testing over 1/8th of the grid showed no significant skill improvement at 500 trees and a degradation at 100 trees. To reduce overfitting with the limited number of historical years, the forest was regularised by limiting tree depth, requiring a minimum number of samples in terminal nodes, randomly sub-sampling predictors at each split, and drawing each tree from a sub-sample of the training data. Candidate configurations were evaluated using nested leave-one-year-out temporal tuning, with the outer held-out year excluded from both hyperparameter selection and final evaluation. The selected reference configuration used a maximum tree depth of 4, a minimum split size of 6, a minimum leaf size of 2, square-root predictor sub-sampling and a bootstrap sample fraction of 0.63. Hyperparameter selection was based on optimizing CRPS relative to the leave-one-year-out climatological reference and average over the individual tests.

Probabilistic forecasts were derived by ranking the predictions of all decision trees at each grid point, and defining tercile thresholds (33rd and 67th percentiles) to represent the probability of onset dates being below-normal, normal, and above-normal. The computational environment included Python (3.7; Van Rossum and Drake 2009), scikit-learn (1.4; Pedregosa et al. 2011), and parallelized computations on an online cluster. Generative artificial intelligence was used to assist with adjusting python code for figure generation (Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8).

The models were trained in a leave-one-out framework over the 24-year training period (1993-2016). For each target year, the model was trained using predictors and onset dates for the remaining 23 years. Hyperparameter selection was performed only within those 23 training years using a nested year-wise cross-validation procedure. The average across year configuration with the highest skill was then used to predict the held-out year. This procedure produced 24 out-of-sample forecasts for each evaluated initialization month. Model skill was evaluated with this nested leave-one-year-out hindcast dataset.

Hindcast skill claims are based on the 1993–2016 leave-one-year-out evaluation, meaning that each target year is temporally out of sample. However, the outer folds share much of their training data and should not be interpreted as independent development and validation datasets. The reported scores instead summarize the complete set of year-wise out-of-sample hindcasts over 1993–2016. The 2023 and 2024 forecasts are presented only as operational case studies.

To test whether the choice of RF was justified, the RF was compared with Bayesian Ridge and multilayer perceptron (MLP) models. All methods used the same leave-one-year-out setup, target variable, grid point mask and verification metrics. The comparison was undertaken using the same predictors. RF uncertainty was represented by the distribution of individual-tree predictions; Bayesian Ridge and MLP predictive distributions were represented using their fitted mean and residual uncertainty. Performance was assessed using CRPS and CRPSS relative to the leave-one-year-out climatological reference. The comparison therefore evaluates the sensitivity of forecast skill to algorithm choice under identical predictor information, rather than identifying a universally optimal model.

3. Results

We assess the forecast performance in three stages. First, we compare the ML methods under a common predictor set and identical leave-one-year-out method. This establishes whether the choice of RF is supported by skill. Second, we evaluate our RF forecast against the raw ECMWF dynamical onset forecast to quantify the added value of statistical post-processing. Finally, we examine the spatial distribution of skill at the country level and illustrate forecast behaviour for recent operational case-study years.

Figure 3 shows positive CRPSS for all three evaluated methods at all four forecast initialisations. The MLP attains the highest Tanzania-average skill for the June, July, and August initialisations (CRPSS of 0.087, 0.095, and 0.122, respectively), while RF has the highest score for the September initialisation (0.114). RF skill is close to that of the other methods in July and August and remains positive from June onwards. Thus, the comparison does not suggest that RF is universally the most skilful method. Its use as the operational reference reflects a combination of competitive skill over Tanzania, stable behaviour across forecast start months, direct ensemble-based probabilistic output, and the possibility of examining the influence of different predictor groups.

We next assess whether RF post-processing provides deterministic added value relative to the raw ECMWF s5.1 onset forecast. This comparison enables the calculation of RMSE for both systems and assess the potential improvement of the RF model over ECMWF forecasts. RMSE was calculated from the ensemble-mean onset forecast and the observed onset date. Figure 4 shows an example of the average difference in RMSE for forecasts initialised in July, averaged over all years. Table 2 shows the RMSE difference between our RF model and ECMWF forecasts for onset date across six countries in Southeast Africa.

Overall, our RF model demonstrates the greatest improvements in Malawi, particularly in August (-7 days) and September (-15 days), and in Mozambique, where the largest reductions occur in July (-7 days) and September (-13 days). Tanzania also shows substantial gains, with errors reduced by 7 days in August and 12 days in September similarly to Zimbabwe that shows errors reduced by –5 in August and –10 in September. Burundi shows modest improvements in July (-3 days) and August (-2 days), but substantial improvements in September (-9 days). By contrast, Rwanda shows little to no improvement, or slight deterioration, with differences of 0–3 days across all months. It is worth noting, however, that Rwanda, and Burundi, have much smaller spatial coverage than the other countries considered.

Seasonally, the strongest improvements are observed in late-season forecasts (August–September), especially in Malawi, Mozambique, and Tanzania. Early-season results (June–July) are, as expected, more mixed due to the longer lead times, with improvements in Malawi and Mozambique, but smaller deteriorations in other countries.

For operational implementation, our RF predictive distribution was represented by 51 tree predictions sampled from the 250-tree forest, matching the 51-member ensemble size of the real-time ECMWF s5.1 seasonal forecast. The historical ECMWF reforecasts available for the 1993–2016 hindcast period contain 25 members; all available members were retained for the hindcast evaluation.

Skill was measured using the Continuous Ranked Probability Score (CRPS) (Hersbach, 2000) and its skill score, CRPSS. For a predictive CDF F and observed value y.

$$\mathit{CRPS}\left(F,y\right)={\int }_{-\infty }^{\infty }{\left[F\left(z\right)-1\left\{z\ge y\right\}\right]}^2\mathit{dz}$$

(3)

With the skill score defined as:

$$\mathit{CRPSS }= 1 - \left({\displaystyle \frac{\mathit{CRP}{S}_{\mathit{forecast}}}{\mathit{CRP}{S}_{\mathit{ref}}}}\right)$$

(4)

Here, the reference forecast used to compute the CRPS_ref is the climatological distribution over the years 1993 to 2016. Positive CRPSS values indicate that the RF model outperforms climatology, while negative indicates deterioration. Our RF model’s performance against climatology was compared to the onset forecast derived from the ECMWF s5.1 dynamical forecasting system performance against the same climatology. Figure 5 shows improvements in CRPSS for forecasts initialised in July, with the largest improvements over the north-eastern part of the domain, especially over Tanzania. Table 3 presents the spatially aggregated differences in Continuous Ranked Probability Skill Score (CRPSS) between our RF model and ECMWF forecasts for onset prediction. Positive values indicate an improvement in forecast skill by our RF model relative to ECMWF.

Table 3 shows the aggregated CRPSS improvements across all countries. A positive RF–ECMWF difference does not necessarily imply positive absolute RF skill relative to climatology, because it can occur when both systems have negative CRPSS but RF has the less negative score. Our RF model generally demonstrates small but consistent improvements, with the largest gains observed in Tanzania (up to +0.15 in August and +0.14 in September) and Burundi (+0.09 in September). Mozambique also see small improvements of +0.02 to +0.05 across all months. Malawi shows only modest increases (+0.01 to +0.04), while Rwanda presents mixed results: small positive improvements in June, July and September, but a notable decline (-0.08) in August. Finally, Zimbabwe shows a similar small improvement for all month (+0.04 to +0.06). Formal significance testing was not undertaken because the number of hindcast years is limited and the data are spatially correlated; the results should therefore be interpreted as indicative differences in predictive performance.

Similarly to the improvements for RMSE, the strongest improvements in CRPSS occur in late-season months (August-September), particularly in Tanzania and Burundi. In contrast, early-season months (June-July) show smaller but consistent skill increases across most countries with exception of Tanzania which for all months see the largest improvements of all the countries.

After running operationally each month as part of the FOCUS-Africa project since early 2023, the model has produced forecasts for two OND seasons, illustrating the range of spatial onset-anomaly patterns represented by the operational product. This is seen in Figure 6, which shows the MSWEP-based forecasts from July of each year alongside the observations. These indicated relatively contrasting large-scale spatial features, with later onset in the southern parts of the domain and earlier onset in the northern regions in 2023, and a reverse pattern in 2024. In both cases, the forecasts reproduce several of the large-scale spatial anomaly patterns seen in the corresponding observations and the estimated actual onset dates from the MSWEP dataset, which displays dipole-like patterns (Figure 7). These two operational cases are illustrative only and, because they cover just two years, do not establish forecast robustness; quantitative performance is evaluated using the 1993–2016 leave-one-year-out hindcasts.

Notably, these forecasts were initialised as early as July 2023 and July 2024, respectively, several months before the OND season. The spatially coherent patterns shown in the July forecasts (Figure 6) and corresponding observations (Figure 7) suggest that large-scale conditions, including ENSO and the IOD, may have contributed to the onset variability observed in these two years. This highlights the potential for useful predictive value several months before the rainy season begins. In fact the dipole-like pattern was evident in all subsequent forecasts issued during both years.

Our RF model’s capability to rank predictors based on their importance provides insights into the factors associated with the onset of the rainy season. Analysis of variable importance reveals that global climate drivers such as ENSO and IOD indices together with local observations make the largest contribution to predictive skill at longer lead times (i.e., June, but also July). As the rainy season approaches (namely in August and September), dynamical model outputs (for all three ECMWF, DWD and CMCC systems) increase their relevance likely due to their ability to capture evolving regional-scale patterns and local atmospheric conditions relevant to onset timing (Figure 8).

4. Discussion and Conclusions

This study demonstrates that statistical post-processing using RF can improve predictions of OND rainy-season onset across much of Southeast Africa relative to raw dynamical forecasts. However, some limitations should be considered when interpreting these results. The 1993–2016 evaluation contains 24 target years, leaving 23 years for model fitting within each leave-one-year-out fold. This period was the longest common hindcast period for which the onset target and predictors from all the dynamical forecasting systems were consistently available. RF complexity was therefore constrained through the nested temporal tuning procedure, and the results should be interpreted as evidence of skill over this specific hindcast period rather than proof of universally robust performance. The many spatial grid points do not constitute additional independent temporal samples, because onset anomalies are spatially correlated. Evaluation over additional independent years, using predictors available in real time, will be required to establish the stability of the system under future climate conditions.

While our RF model performed well across most of the study area, certain local areas, such as western Zimbabwe, exhibited lower skill. This may be due to local climatic complexity and/or data sparsity. Future improvements could include incorporating additional local predictors such as soil moisture, vegetation indices, while the inclusion of sub-seasonal forecasts represents another promising direction for future research. Further increasing the resolution of observational datasets could also help make the model more relevant to the end-users. Furthermore, the definition of an OND onset season used is tailored to Tanzania, reflecting its bimodal and unimodal rainfall regimes. This may partly explain the larger skill improvements observed there compared to other countries. This underscores the need for continued stakeholder engagement and collaboration with local meteorological agencies to refine model definitions and enhance operational relevance.

Improved understanding of the physical mechanisms linking ENSO and IOD to regional rainfall variability could further enhance model performance. This could include exploring additional lead or lag times for some predictors over certain areas, and may also help anticipate how future climate change-driven shifts in these large-scale drivers might affect seasonal rainfall patterns.

In this study, RF is not assumed a priori to be the most skilful machine learning method. It was selected following a comparison with Bayesian Ridge and multilayer perceptron models under the same leave-one-year-out hindcast framework and common predictors. Although MLP achieved higher skill in some start months, RF was retained as the main model because it exhibited more consistent performance across all start months, it provides an ensemble-based probabilistic forecast directly and its behaviour can be easily analysed through its predictor-importance diagnostics.

Our RF model’s ability to integrate multiple predictors—including observational indicators and outputs from dynamical forecast models—led to an improvement over ECMWF’s s5.1 dynamical forecasts across much of the study region, particularly for August-September forecasts, confirming its potential as a powerful tool for operational forecasting. Utilizing forecasts from multiple dynamical models provided complementary information and may have allowed the RF approach to leverage the strengths from different systems across different areas of the domain, contributing to improved predictive performance.

Our RF model achieved measurable improvements in prediction accuracy, with notable regional differences that may partly reflect the influence of large-scale climate drivers such as ENSO and IOD. In particular, our RF model captured contrasting north-south onset dipole patterns in 2023 and 2024 that were associated with differing large-scale climate conditions, demonstrating its ability to translate large-scale variability into actionable local insights. With continued refinement, such a hybrid modelling approach can help mitigate the impacts of climate variability by providing timely and reliable information to decision-makers in sectors such as agriculture, water management, and disaster preparedness.

To quantify the improvement, our RF model reduced RMSE by 2 to 15 days compared with ECMWF s5.1 forecasts, with the largest gains observed for the August-September forecast initialisations—the key decision-making window ahead of the OND rains. Probabilistic skill, measured via CRPSS, improved by up to 0.15 over the same reference system. These gains highlight the operational relevance of the RF approach, especially in leveraging large-scale climate drivers at longer lead times and incorporating dynamical forecast strengths closer to the rainy season, thereby providing a more reliable and actionable forecast tool for Southeast Africa.

While this study demonstrates potential for application beyond the study domain, further validation is necessary to adapt the RF methodology effectively across different climatic regions. Importantly, the FOCUS-Africa project has demonstrated that strategic and sustained collaboration between international research groups, meteorological agencies, and local stakeholders is essential to maximize the practical utility and impact of advanced forecasting tools like RF. Such collaborative efforts are crucial for enhancing regional resilience against climate variability and supporting sustainable socio-economic development across Southeast Africa. As a direct outcome of the project, TMA is now using an operational version of the onset prediction system, along with forecasts of cessation, rainfall and temperature. These forecasts are delivered through a tailored version of the Teal tool (https://tealtool.earth) developed by the World Energy & Meteorology Council (WEMC). Figure 9 illustrates the delivery of such seasonal forecasts, showing onset predictions for different regions across the domains for August 2024.