SPEAR Model and Experiments.

The main conclusions of this study are obtained from the SPEAR (33, 42). SPEAR represents the newest modeling system for seasonal to multidecadal prediction which incorporates new model development components that have occurred in the last decade at the National Oceanic and Atmospheric Administration (NOAA) GFDL. These include: a new dynamical core (43), revised atmospheric physics (44), a new sea-ice and ocean model (45), and an enhanced land model (46). The SPEAR atmospheric model uses 33 levels in the vertical and is run at different atmospheric-land horizontal resolutions: 0.5° (SPEAR_MED) and 1° (SPEAR_LO) in this paper. The intermediate 0.5° configuration, SPEAR_MED, is a compromise between the possibility to run a large ensemble of simulations with available computation resources and retaining enough horizontal resolution to study regional climate and extremes. It is worth noting that the SPEAR_MED large ensemble features a horizontal grid spacing (0.5°) that is finer than those used in most of the previously used large ensembles [with the exception of FLOR (31)], thus making these GFDL ensembles a unique and unprecedented tool to study extremes and regional climate.

We used four different numerical experiments: 1) a long-term control simulation (CTRL) to evaluate unforced natural variability; 2) an ensemble driven by natural forcing only (NATURAL) to provide a baseline with only natural forcing (i.e., volcanic eruptions and solar cycles); 3) an ensemble driven by observed natural and anthropogenic forcing up to 2014 (HIST) and then according to the intermediate (≈+3 K of global warming by the end of the 21st century) SSP2-4.5 developed for the CMIP6 (35, 47); and 4) an ensemble driven by observed natural and anthropogenic forcing up to 2014 (HIST) and then according to the CMIP6 high-emission, fossil fuel-dominated (≈+5 K of global warming by the end of the 21st century) SSP5-8.5.

The 3,000-y CTRL simulation is driven by CO2 forcing kept constant at 1850 levels. Climate drifts associated with this long-term integration are estimated to be very small and statistically insignificant for the winter SA rainfall. The 30 members of the NATURAL ensemble are driven by the same observed natural forcing (i.e., solar and volcanic) until year 2014, and by only solar forcing (quasi-11-y cycle) after 2014, with the anthropogenic forcing held fixed at the 1921 level. In the HIST + SSP5-8.5 (HIST + SSP2-4.5) ensemble, each member is driven by observed natural and anthropogenic forcing (greenhouse gases, anthropogenic aerosols, and ozone) up to year 2014 and by the SSP5-8.5 (SSP2-4.5) forcing afterward. More details about how the SPEAR large ensemble is designed can be found in Delworth et al. (2020) (33).

Model Evaluation.

In addition to the model’s ability to reproduce the wintertime southern African climatology (Fig. 1 C–E), the performance of SPEAR_MED in simulating wintertime rainfall variability and historical trends (1951–2017) over SSA was evaluated against three different observational land rainfall datasets: the Global Precipitation Climatology Center (GPCC) dataset (48), version 7; the Climate Research Unit high-resolution grids of monthly rainfall at the University of East Anglia (49), version 3.24; and the University of Delaware (UDEL) precipitation dataset, version 5 (http://climate.geog.udel.edu/∼climate/), all at 0.5° resolution. The choice of these three gridded observed datasets, in place of scattered measurements from the South African Weather Service meteorological stations, was dictated by the need to be able to compare models with observations, as done in previous studies (5). The values of these three precipitation datasets for SSA were obtained from a limited number of stations and different interpolation algorithms. As a consequence, they can feature, locally, considerable differences (e.g., Fig. 1A and SI Appendix, Fig. S1). However, differences in area-averaged metrics like, e.g., the WRI, are minimal (Fig. 1B), thus making our results independent from the choice of the single precipitation dataset.

In order to have a realistic representation of the width of the distribution of rainfall anomalies, it is key that SPEAR_MED reproduces the interannual, multiannual, and decadal natural variability of the SSA winter rainfall. To check this, we worked out the SD of the detrended full, 3-y, and 10-y low-pass-filtered WRI from the three observational datasets and the SPEAR_MED ensemble members over the common period 1921–2010 (SI Appendix, Fig. S3). The SD of the observations was between 5 mm⋅month−1 (CRU) and 6 mm⋅month−1 (GPCC and UDEL) for the 3-y low-pass-filtered WRI. The SD values from the model ranged from 4 to 6.3 mm⋅month−1. The observed values are, therefore, within the range from the model, suggesting that the model has the ability to properly estimate the magnitude of 3-y lasting droughts. Similarly, a good agreement between SPEAR_MED and observations exists for the SDs calculated from the unfiltered WRI time series (interannual variability) and from 10-y low-pass-filtered WRI (decadal and longer variability), too.

The effect of internal natural variability is large for SSA winter rainfall (27⇓–29); thus, it is not appropriate to compare observed AMJJAS rainfall trends directly with the ensemble mean or with each single ensemble member, which may show contrasting signs (SI Appendix, Fig. S2). Instead, we evaluated if SPEAR_MED’s historical trends of AMJJAS rainfall are consistent with observations over SSA. To do so, we computed rainfall trends over the last 67 y (1951–2017) in GPCC, CRU, and UDEL and compared them with individual members of the HIST + SSP5-8.5 ensembles over the same time period.

If the observed trend at one grid point was within the range of those simulated by the 30 HIST ensemble members, then we said that the model is consistent with observations in that grid box. We found that SPEAR_MED is consistent with observations over most of southern Africa (SI Appendix, Fig. S13).

Additional Large Ensembles.

To assess the robustness and model dependence of our results, we analyzed five additional large ensembles (SI Appendix, Table S1): 1) the SPEAR_LO ensemble (33); 2) the GFDL FLOR model, at 0.5° land/atmosphere resolution; 3) the flux-adjusted FLOR (FLOR_FA) large ensembles, obtained imposing temperature and salinity flux adjustments at the ocean surface to FLOR (50) (both with a land-atmospheric horizontal resolution of 0.5°); 4) the CESM-LENS (30), with land-atmospheric horizontal resolution of approximately 1°; and 5) the MPI-GE (26), with land-atmospheric horizontal resolution of 1.8°. These additional large ensembles are available with various CMIP5 scenarios and are documented in SI Appendix, Table S1. An evaluation of the wintertime climatology over SSA showed that these models all underestimate AMJJAS mean rainfall (Fig. 1 C–E and SI Appendix, Fig. S14 and Table S2). With the exception of SPEAR_LO, these models also underestimate the SD of the full 3-y and 10-y low-pass-filtered WRI (SI Appendix, Fig. S3). Critically, this means that they also underestimate the width of the probability distribution of the 3-y AMJJAS rainfall anomalies (SI Appendix, Fig. S15). In particular, CESM-LENS and FLOR_FA have SDs that are 50% and 40% smaller, respectively, suggesting that they are poor tools for risk analysis over SSA. As they substantially underestimate the probability of occurrence of event_1517, to quantify changes in risk in a manner that implicitly accounts for model biases, we used a 3-y WRI anomaly corresponding to the first percentile, which is the percentile to which −11.5 mm/month corresponds to in observations and SPEAR_MED.

WRI.

In this study, we focused on the regional-scale drought of the Western Cape. We thus used the annual time series of the WRI (29) to monitor interannual variability and monthly rainfall anomalies. To define the WRI, we first selected the grid points where at least 65% of the total annual rainfall occurred from April to September (Fig. 1 C–E and SI Appendix, Fig. S14). Then, we took the areal mean of the extended winter (i.e., April–September) rainfall over the irregular region defined above (Fig. 1 C–E and SI Appendix, Fig. S14). The WRI is, thus, the area-averaged rainfall over the portion of SSA that experiences a dry summer and a wet winter—that is, a Mediterranean rainfall regime. This area encompasses the region of intensely irrigated agriculture surrounding the metropolitan area of Cape Town, as well as the water basins of the Breede and Berg Rivers, where dams supplying water to Cape Town are located.

Detectability of the Mean Rainfall Change.

To determine where and when the decadal changes in AMJJAS rainfall start being caused by external forcing, and not by multidecadal variability, we applied a Monte Carlo approach to the long CTRL run: At each grid box, we randomly chose a 10-y period and a nonoverlapping 50-y period (to mimic 1921–1970). Then, we computed the time mean difference between the 10-y and 50-y time series. This difference is solely associated with internal natural variability of the climate system. This process was repeated 30 times (to mimic the 30-member ensemble); we then took the ensemble mean of these differences. The whole process was then repeated 10,000 times to create an empirical probability distribution of these ensemble mean differences, which was used to assess the detectability of decadal changes in rainfall. Anomalies outside the range of the distribution were attributed to external forcing and considered detectable against internal climate variability (Fig. 4 and SI Appendix, Figs. S11 and S12).

Estimation of the Probability Distribution.

We derived a probability distribution of the 3-y mean WRI anomalies due to natural variability alone from the long CTRL run. We randomly selected a 50-y and 3-y sequence (nonoverlapping) and then calculated the anomaly of the 3-y period relative to the 50-y climatology. This choice mimics the 2015–2017 mean minus the 1921–1970 mean. We repeated this process N times (N = 10,000) to form a distribution of the 3-y WRI anomalies (Fig. 2 A–D). The probability of occurrence of experiencing a 3-y WRI anomaly equal to or less than the 2015–2017 anomaly—as per the gridded datasets—is about 1% in CTRL and 0.7% from HIST, taking the average of decadal probabilities over 1921–1970, respectively (Fig. 2E). Similarly, we estimated the distribution of the 4-y WRI anomaly. The probability of occurrence of a WRI anomaly of the same intensity, but of one additional year of duration, is 0.4% and 0.2% from the CTRL and HIST, respectively.

To evaluate the decadal change in the probability of occurrence of a 3-y WRI anomaly equal to or worse than that of 2015–17, we empirically estimated a decadal-varying probability distribution using the HIST and SSP5-8.5 (SSP2-4.5) experiments. The probability distribution was estimated for a 20-y time window, so that, for example, that referred to as 2010 is built from all years from 2001 to 2020. This choice was motivated by the need to have a time period not too wide in order to assume the stationarity of the probability distribution, but, at the same time, a number of instances large enough to allow for sufficiently accurate estimates of probabilities of rare events (e.g., 100-y return time). In a 20-y time window, there are 18 different 3-y WRI anomalies (relative to the climatological reference period 1921–1970). This leads to 18 × 30 = 540 different values when considering all of the 30 ensemble members, from which we empirically built the decadal probability distribution. Once we have the decadal probability distribution, we can estimate the probability of occurrence, for each bidecadal period, of 3-y WRI anomaly equal to or less than that observed in 2015–2017 (−11.5 mm⋅month−1, obtained averaging GPCC, CRU, and UDEL) for any random 3-y segment within the 20-y time window. The 95% CIs in these probabilities were estimated by applying bootstrap-with-replacement resampling 10,000 times. The same methodology was applied to estimate the probability of occurrence of 4-y droughts.

We quantified the uncertainty in the estimate of the decadal probability of occurrence, derived from only 540 different 3-y rainfall anomaly values, as follows: We took the long 3,000-y CTRL and randomly selected 50-y and 3-y nonoverlapping periods and estimated the difference. We repeated this step N times (with N = 10,000) to obtain a large population sample of N 3-y anomalies, from which the probability of the event_1517 was estimated to be ≈1%. From this large sample, we then randomly drew M realizations (with replacement), with M≤N, and estimated the probability of occurrence. For each value of M, we repeated the last step 10,000 times and obtained 10,000 different probability estimates, which allowed us to estimate the 95% CI (SI Appendix, Fig. S4B). As expected, the CI decreased with M up to approximately [0.9%, 1.2%] for M = 10,000. For values of M less than 300, the uncertainty is so large that it is impossible to have any sensible estimate of the probability of the event. For M = 540, the CI is approximately [0.5%,1.7%], which we can consider sufficiently accurate for our purposes.

Joint Probability Distribution of Drought Intensity and Duration.

The probability distribution of a drought in the Cape Town’s Mediterranean area as a function of duration and intensity was estimated from the historical and projected AMJJAS WRI anomaly time series. The focus in this paper is on severe droughts; therefore, we selected, for each time series, all contiguous years for which the WRI anomaly was below −0.75 SD (≈ −6 mm⋅month−1). With this choice, we excluded years that were moderately and very moderately dry. For each of these segments, we worked out the mean WRI anomaly by averaging the annual WRI anomaly values over the whole segment. We chose a 2-mm⋅month−1× 1 y bin (Fig. 3) to work out the percentage of the droughts within each bin. The analysis was performed for the 1921–1970 time period and for the periods 2011–2040, 2041–2070, and 2071–2100. To evaluate if the probability differences relative to 1921–1970 are attributable to anthropogenic climate change, we applied the same method to the 3,000-y CTRL. We randomly selected a 50-y and a 30-y nonoverlapping time span and computed the number of droughts for each duration-drought intensity bin. We repeated this 30 times to mimic the 30-member ensemble and so worked out the probability differences between the 50-y and 30-y periods. The whole process was then repeated 10,000 times to create an empirical probability distribution of the probability differences for each bin: Anomalies outside the range of the distribution were attributed to external forcing and considered detectable against internal climate variability.