Projected changes in precipitation characteristics over the Drakensberg Mountain Range

This study examines the potential impacts of climate change on the characteristics of precipitation over the Drakensberg Mountain Range (DMR) at different global warming levels (GWLs: 1.5, 2.0, 2.5 and 3.0°C) under the Representative Concentration Pathway 8.5 (RCP8.5) scenario, using dynamical and statistical downscaled datasets. The dynamical datasets consist of 19 multi‐model simulations datasets from the Coordinated Regional Climate Downscaling Experiment (CORDEX), whereas the statistical downscaled datasets comprise 19 multi‐model simulations from the National Aeronautics and Space Administration (NASA) Earth Exchange (NEX) Global Daily Downscaled Projections (NEX‐GDDP, hereafter NEX). The capacity of the CORDEX and NEX datasets to represent past characteristics of extreme precipitation over the DMR was evaluated against eight observation datasets. The precipitation characteristics were represented by eight precipitation indices. Both CORDEX and NEX realistically capture the characteristics of extreme precipitation over the Drakensberg and, in most cases, their biases lie within the observation uncertainty. However, NEX performs better than CORDEX in reproducing most of the precipitation characteristics, except in simulating the threshold of extreme rainfall. The ensemble means of CORDEX and NEX agree on a future increase in the intensity of normal precipitation, in the frequency and intensity of extreme precipitation, as well as an increase in widespread extreme events, with a decrease in the number of precipitation days and continuous wet days. However, they disagree on the projected changes of annual precipitation, for which CORDEX projects an increase over most parts of the DMR, whereas NEX indicates a decrease. The self‐organizing‐map analysis, which reveals diversity in the projection patterns hidden in the ensemble means, shows the most probable combinations of projected changes in the annual precipitation and extreme precipitation events (in terms of intensity and frequency): (a) increase in both annual precipitation and extreme precipitation events; (b) decrease in both annual precipitation and extreme precipitation events; (c) decrease in annual precipitation but increase in extreme precipitation events. The results of this study can thus provide a basis for developing climate change adaptation and mitigating strategies over the DMR.


| INTRODUCTION
The Drakensberg Mountain Range (DMR) is one of the most valuable natural resources in Southern Africa. Its natural beauty makes it one of the major tourist attractions of the region. The DMR features pristine steep-sided river valleys, rocky gorges, numerous caves and rock shelters (which contain about 665 rock art sites). It provides a natural habitat for many plants and animal species. For example, it is home to endemic and threatened species like the yellow-breasted pipit. Due to this rich natural diversity, DMR was designated a UNESCO World Heritage Site in 2000 (UNESCO, 2000). Furthermore, the DMR is the source of rivers that support socio-economic activities, like agriculture, mining and hydro-electric power generation (Tyson et al., 1976;Nel and Sumner, 2006). For example, the Vaal River, whose source lies in the upper catchments of the Drakensberg, is the largest tributary of the Orange River-a key water source for industrial activities in the dry inland of South Africa and Namibia. The DMR also supports two large inter-basin water transfer schemes: the Tugela-Vaal water transfer scheme and the Lesotho Highlands Water Project. Hence, rainfall variability over the DMR poses a threat to the environment and water resource management in Lesotho, South Africa and Namibia, because a substantial decrease in the annual rainfall can lead to hydrological droughts in all these countries. Conversely, an extreme rainfall event over the DMR can cause floods, landslides and soil erosion, as well as the destruction of properties and loss of lives in the communities around the DMR (Hyden et al., 2000). For instance, an extreme rainfall event in 1987 led to a flood that caused more than 300 fatalities and damaged infrastructure worth a billion Rand (about 85 million U.S. Dollars) in the KwaZulu-Natal Province of South Africa (Tennant and van Heerden, 1994). The recent KwaZulu-Natal and Eastern Cape floods and landslides (11-13 April 2022) illustrate the devastating impacts floods can have on life and property. The floods, which battered towns in the region, left 4,000 houses destroyed, 40,000 people displaced, 48 people missing and at least 443 dead (ECHO, 2022). Meanwhile, there is concern that the ongoing global warming may enhance the intensity of rainfall variability over the DMR (Mason et al., 1999;Karl and Trenberth, 2003;Kruger, 2006;Fatti and Patel, 2013;Abiodun et al., 2017;Pohl et al., 2017). Hence, to reduce the socioeconomic impacts of rainfall variability in this area, there is a need to study the potential impacts of climate changes on extreme precipitation indices over this mountain range.
Mountain environments are very sensitive to climate variability and change because of their complex topography (Thompson, 2000;Beniston, 2003;Diaz et al., 2003). The highly variable topography of mountains with steep slopes causes climate variables (e.g., temperature and precipitation) to change rapidly and systematically with height over relatively short horizontal distances (Whiteman, 2000). This is also the case with vegetation and hydrology variables (Becker and Bugmann, 1997). Hence, mountains exhibit high biodiversity, often with sharp transitions (Beniston, 2006) in vegetation sequences and rapid changes from vegetation and soil to snow and ice. However, the high relief and sharp gradients also make mountain ecosystems very vulnerable to slight changes in temperature and extreme precipitation events. Diaz et al. (2003), for instance, showed that a relatively small perturbation in global processes can cascade down to produce large changes in most or all the myriad interdependent mountain systems, from the hydrological cycles to the complex fauna and flora, and down to the people who depend on those resources. Beniston (2003) also showed that mountainous areas that experience more favourable temperature and moisture conditions, as a result of global warming, may become breeding grounds for mosquitoes, increasing the spread of malaria. Changing climates may furthermore alter the seasonal patterns of tourism over mountains (e.g., skiing in winter), and thus the environmental pressure associated with different forms of tourism today. Consequently, given that mountains differ considerably from one region to another, there is a need to investigate the impact of climate change on individual mountains. While many studies have discussed the impacts of global warming on various parts of southern Africa (Hulme et al., 2001;Kruger and Shongwe, 2004;Warburton et al., 2005;New et al., 2006;Kusangaya et al., 2014;Pinto et al., 2016;Mbokodo et al., 2020;Barimalala et al., 2021), there is a dearth of information on how global warming could alter the characteristics of precipitation over the Drakensberg, especially at various global warming levels (GWLs). Such information would be very useful for policymakers in all countries whose socio-economic activities depend on DMR resources.
General circulation models (GCMs) are the principal tool for projecting future climate. However, due to their coarse horizontal resolution, they are unable to resolve subgrid (or local) scale meteorological processes, which affect regional climates and potentially have strong feedback on the regional climate (Wilby and Wigley, 1997). Hence, GCM projections are inadequate for studying the impacts of climate change over complex topography such as the DMR, and how these local scale processes will be affected in the projected climate. To overcome this shortcoming, GCM simulations are usually downscaled to the local scale to incorporate the impacts of small-scale features. There are two major approaches for downscaling GCM outputs: statistical and dynamical downscaling approaches. Statistical downscaling is based on establishing a linear or nonlinear relationship between subgridscale parameters and grid-scale (coarser) resolution predictor variables (Wilby and Wigley, 1997) while dynamical downscaling is based on nesting a limited area model (LAM) in a GCM, such that the GCM fields provide initial and lateral boundary conditions for the LAM simulation Mearns, 1991, 1999). Both approaches are associated with some uncertainties. For example, inadequate representation of local scale features, such as the topography in a GCM, might lead to the introduction of errors in the simulated meteorological fields (Giorgi and Mearns, 1999). In addition, nesting LAMs results in a sharp transition between the coarse and high resolutions and introduces spurious numerical noise (Giorgi and Mearns, 1999). This can distort wave propagation and reflection properties (Qian et al., 1999), thereby resulting in errors in the simulated fields. Conversely, statistical downscaling does not account for the feedback from local scale features (topography, land use and land cover change), which can either amplify or attenuate the projected climate change signal (Hewitson and Crane, 2006;Thrasher et al., 2012). Furthermore, the statistical downscaling approach assumes that the relationship between the large-scale and the local-scale climate in the present climate will remain the same (stationary) in a projected future climate (Wilby and Wigley, 1997;Hewitson and Crane, 2006). The validity of such an assumption cannot be fully assessed. Given the uncertainty associated with both downscaling approaches, a comparison of projections from both methods will help improve the robustness of climate projections at the local scale. Most climate change projection studies over South Africa (i.e., Tadross et al., 2005;Engelbrecht et al., 2009;Pinto et al., 2016;Maure et al., 2018) have used dynamical rather than statistical downscaling. However, to arrive at a more robust climate change projection, there is a need to compare the results of both approaches. Recently, Landman et al. (2018) and Abiodun et al. (2020) employed both statistical and dynamical approaches to study the impacts of climate change over southern Africa, but the focus of their studies was not over the DMR.
Hence, the aim of this study is to examine the impacts of climate change on the characteristics of precipitation over the DMR, using both dynamically and statistically downscaled datasets. The impact of climate change is considered at various GWLs under RCP8.5. Section 2 provides the sources and characteristics of the datasets, definitions of the precipitation indices and methods used in analysing the datasets. Section 3 presents and discusses the results of the study, while Section 4 gives the concluding remarks.

| Study area
This study focuses on the DMR (26 -32 E; 28 -32 S; Figure 1). With the highest peak at about 3,475 m above mean sea level, the DMR is the highest mountain range in Southern Africa (Figure 1). This trans-frontier mountain range straddles the border between South Africa and Lesotho. This region also receives its highest rainfall in austral summer and the lowest in austral winter. Thunderstorms and orographically induced storms are major sources of rainfall over the DMR (Tyson et al., 1976). Cold fronts that move across southern Africa from a west-northwesterly to east-southeasterly direction also bring widespread rainfall and occasionally snow over the Drakensberg (Tyson et al., 1976). Summer rainfall variability over the region at sub-seasonal scales is influenced by systems such as Tropical Temperate Troughs (Hart et al., 2013) and mesoscale convective complexes (Blamey and Reason, 2013). Inter-annual rainfall variability is influenced by El Nino Southern Oscillation (ENSO; Nel and Sumner, 2008;Dieppois et al., 2016). Decadal and inter-decadal summer rainfall variabilities have been found to be associated with changes in the meridional circulation between the South Atlantic and Indian Oceans, which may relate to variations in the subtropical ridge of the first zonal standing wave over the Southern Hemisphere. According to Malherbe et al. (2014Malherbe et al. ( , 2016, this could be consistent with atmospheric anomalies related to the Southern Annular Mode (SAM) during the late austral summer (Dieppois et al., 2016).

| Data
Eight observation datasets and two simulation datasets were analysed in this study. Information about the eight observation datasets is provided in Table 1. Some studies (e.g., Sylla et al., 2013;Abiodun et al., 2017Abiodun et al., , 2020Kumi and Abiodun, 2018) have analysed these datasets and documented the usefulness and reliability of the datasets Note: The period of the data analysed in each dataset is indicated.
over Africa. In this study, the datasets were used to evaluate the performance of the simulation datasets in reproducing the characteristics of rainfall over the Drakensberg. The simulation datasets used were the Coordinated Regional Climate Downscaling Experiment (CORDEX) and the National Aeronautics and Space Administration (NASA) Earth Exchange Global Daily Downscaled Projections (NEX-GDDP, hereafter NEX; Thrasher and Nemani, 2015). NEX has a spatial resolution of 0.25 , and its temporal coverage spans the period 1950-2100. The NEX simulations dataset consists of statistically downscaled climate scenarios for the globe, obtained from General Circulation Models (GCM) runs of the Fifth phase of the Coupled Model Intercomparison Project (CMIP5; Taylor et al., 2012). It includes 42 climate projections from 21 CMIP5 GCMs and two Representative Concentration pathways (RCPs; Van Vuuren et al., 2011) scenarios (RCP4.5 and RCP8.5). Only the RCP8.5 is analysed for this study (see list of models in Table 2). We used the high emission scenario (i.e., RCP8.5) dataset because it has the largest number of simulation ensemble members in both datasets. For most models, the RCP4.5 simulation does not attain 2.5 or 3.0 C global warning levels. Using the fewer available RCP4.5 simulations will compromise the comparison at these two warming levels. The NEX climate change projections are bias-corrected using the Bias-Correction Spatial Disaggregation (BCSD) Method (Wood et al., 2002(Wood et al., , 2004. The BCSD algorithm utilizes the spatial detail provided by observationally derived datasets to interpolate GCM outputs to higherresolution grids. This makes the NEX data suitable for T A B L E 2 The GCMs simulations, corresponding 30-year period for various GWLs (1.5, 2, 2.5 and 3.0 C) under the RCP8.5 scenario, and the models used for downscaling the GCM simulations used in the study.

GCMs
Period of GWLs Downscaling models 1.5 C 2 C 2.5 C 3 C 1995-2024-2038-20492031-2060 studying climate over complex terrain regions like the Drakensberg. Some of the studies that have used the NEX data include Abiodun et al. (2020); Bao and Wen (2017); Chen et al. (2017) and Yu et al. (2018). The COR-DEX data used for this study were obtained from downscaled simulations of the CMIP5 experiment. Specifically, the CORDEX Africa Regional Climate Model (RCM)-GCM matrix data described in Nikulin et al.

| Definition of timing for GWLs
This study adopted the definition of timings for GWLs used by . Specifically, the timing of the GWL is the first year of a 30-year period, when the global temperature is above the desired warming level (i.e., 1.5, 2.0, 2.5 and 3.0 ) relative to pre-industrial levels (which are defined as 1861-1890). The timing of the GWL can either be defined using only GCMs (Nikulin et al., 2018) or using both GCMs and observation (Vautard et al., 2014;Dosio and Fischer, 2018). Both approaches consist of choosing a baseline period (i.e., an averaged window period of 15, 20 or 30 years, which is usually from the pre-industrial period) and computing the temperature difference between the baseline period (pre-industrial) and the period when the specific GWL of interest (1.5 C or higher) is reached (Vautard et al., 2014;Dosio and Fischer, 2018). However, they differ in that the approach, which uses both observation and GCM data, computes the GWL by calculating the observed global temperature rise between the pre-industrial period and the present period, before adding this to the GCM projected future warming relative to the present. This approach may result in less spread in GWLs across the models, as all GCMs are brought to the same level of warming relative to the present climate. This is because the approach assumes equal climate sensitivities across all GCMs from the pre-industrial periods to the present. It may also be affected by observational uncertainty and artificially reduced/enhanced GCM climate sensitivities . In contrast, the GCM approach which was used in this study and in Nikulin et al. (2018) determines GWLs from GCM experiments that are run from the pre-industrial period through to the 21st century. This provides a spread of GWLs across models (see Table 2) as a result of different climate sensitivities across the GCMs. This approach defines the pre-industrial period as used in all the historical simulations of the CMIP5 models. For each downscaling approach used in this study, we have used the same GWL timing as defined by the corresponding driving GCM to extract a 30-year period for our analysis.

| Rainfall index and climate change computations
The indices used in this study to analyse the characteristics of extreme rainfall events over the Drakensberg are similar to those used by Abiodun et al. (2020). Only T A B L E 3 Description of extreme rainfall indexes used in the study. extreme rainfall indices that could help to quantify the frequency, duration and intensity of extreme events were analysed. The indices are listed in Table 3. For each downscaled simulation (CORDEX and NEX), we computed each index during the control period  as well as during the 30-year period when the simulation (GCM) reached a given GWL. The GWL varies from one GCM to another. We chose the period  as the control period because we want our results to be comparable with other CORDEX Africa studies (Maure et al., 2018;Abiodun et al., 2019) and climate impact studies (e.g., Sakalli et al., 2017) that used the period as their control period. To determine the impact of climate change at different GWLs, we subtracted the simulation of the reference climate from that of the GWL. This was done for each grid point within our study domain ( Figure 1). Following Abiodun et al. (2020) and , the significance of the change was assessed by using two conditions. The first condition used a t-test to check that the climate change at a given grid point was statistically significant for at least 80% of the simulations (with respect to the inter-annual variability of the reference climate [present-day]). The second condition checked that at least 80% of the simulations agreed on the sign of the change at the grid point. The climate change signal at every grid point was considered robust when both conditions are met.

| Testing the robustness of projected changes in spatial patterns
To test the robustness of the projected changes in spatial patterns of the rainfall indices, Self-Organizing Maps (SOM) were used to classify the simulations into groups. A SOM is a type of Artificial Neural Network (ANN) that is used for dimensionality reduction. It takes in high-dimensional data as input and reduces it to a lowdimensional (typically two-dimensional) input space of training samples called maps (nodes). Each node represents a range of spatial patterns within the continuum described by the original data space. A SOM is trained using unsupervised learning and can be used to identify features in data (Kohonen, 1990). Several studies have demostrated the use of SOM in grouping climate variables (e.g. Abiodun et al., 2020;Abatan et al., 2023). In this study, the SOM was used to classify the projected spatial patterns of four extreme rainfall indices, the rainfall total (RTOT), the number of wet days (WDAYS), the extreme rainfall threshold (R97.5p), and the rainfall frequency (R97.5pFREQ) at four GWLs into 12 groups (i.e., a 4 × 3 nodes). We tested the sensitivity of our SOM classification to the choice of 12 groups (4 × 3), 9 groups (3 × 3) and 16 groups (4 × 4) and found that the general patterns of the classification do not change much. However, the 12 groups classification gives a better distinction of the pattern than the 9 groups, whereas the 16 groups classification produces a node with no frequency, indicating too many numbers of groups. For the SOM analysis, the NEX dataset was first re-gridded onto the CORDEX grid. Then, both CORDEX and NEX datasets were combined to form a single dataset, which was used in the SOM analysis. However, the SOM analysis was performed for individual precipitation index, and the contributions of the CORDEX and NEX simulations to the SOM nodes were determined at the four GWLs.
The analysis helped to reveal the diversity in the projection patterns that are usually hidden in ensemble mean results.

| Evaluation of CORDEX and NEX datasets
Here, we evaluate the capability of the CORDEX and NEX simulations to reproduce the characteristics of the rainfall indices (Table 3) over the DMR by comparing the simulation results with the observation results. The evaluation focuses on how well the datasets capture the rainfall intensity and frequency distribution, as well as the mean value and spatial distribution of the rainfall indices over the Drakensberg. However, to put the model evaluation into the right perspective, we first present the differences (i.e., uncertainties) in the results of the eight observation datasets (ARC, CHIRPS, PER-SIANN, TAMSAT, AgCFSR, AgMERRA, WFDEI-CRU and WFDEI-GPCC) before discussing the evaluation of the model. There is good agreement among the eight observation datasets on the precipitation intensity-frequency curves over the DMR, but there are some discrepancies as well (Figure 2). While all the observation curves agree on a decrease in the rainfall frequency with an increase in the rainfall intensity, the rate of the decrease differs among the datasets. The decrease is fastest in PER-SIANN and TAMSAT (which lie at the lower edge of the observation spread) and slowest in AgCFSR (which lies at the upper edge of the observation spread). This implies that, over the DMR, AgCFSR reports the most frequent heavy rainfall, while PERSIANN and TAMSAT report the least frequent heavy rainfall. The tail of the AgCFSR curve extends beyond 120 mm day −1 , whereas that of TAMSAT truncates below 80 mm day −1 . Other datasets cluster together between these two extremes.
The differences between AgCFSR and TAMSAT may be attributed to their data sources. TAMSAT is a satellite product, while AgCFSR is a reanalysis product (Table 1).
It is difficult to determine which of these two sources is better. This is because cloud cover and orbiting of satellites could make the satellite product miss a heavy F I G U R E 2 Precipitation intensityfrequency over the Drakensberg Mountain Range as depicted by (a) observation, (b) Coordinated Regional Climate Downscaling Experiment (CORDEX) and (c) National Aeronautics and Space Administration (NASA) Earth Exchange (NEX) [Colour figure can be viewed at wileyonlinelibrary.com] rainfall event or underestimate the intensity of the event. Conversely, the interpolation and extrapolation of the data in reanalysis could also make the reanalysis product overestimate the intensity of heavy rainfall. Several studies have reported a similar discrepancy among observation products over Africa (e.g., Sylla et al., 2013;Abiodun et al., 2020). Abiodun et al. (2020) even found that the outlier datasets vary over southern African cities. For example, they showed that, while PERSIANN and CHIRPS were outliers over Cape Town, TAMSAT and PERSIANN were outliers over Lilongwe, and AgCFSR over Johannesburg. These results make it difficult to isolate outlier datasets from the model evaluation. Hence, there is a need for more studies on how to reduce disparities among the observation datasets used for model evaluation.
Both simulation datasets (CORDEX and NEX) give a realistic representation of the precipitation-intensityfrequency curves, but among the two, NEX provides a better representation (Figure 2). With CORDEX (Figure 2b), about half of the simulated curves fall within the observation spread, while the other half produce longer tails than those of the AgCFSR curve (i.e., they simulate more frequent heavy rainfall events). With NEX (Figure 2c), however, only four of the simulated curves (i.e., 25%) fall outside the observation spread (i.e., have a longer tail than that of the AgCFSR curve). For the mean values of the precipitation indices over the DMR (Figure 3), the spread of the NEX simulations is closer to the observations than those of the CORDEX simulations. For instance, all the NEX simulations fall within ±50% threshold of the observed values in six indices (i.e., RTOT, WDAYS, SDII, R95.7p, R95.7pTOT and Rx5day), and more than half of the models fall within the threshold in the three indices (CWD, R20mm and WERE). In contrast, all the CORDEX simulations only fall within this threshold in one index (i.e., SDII) and more than half of them fall within this threshold in four indices (WDAYS, CWD, R95.7 and WERE), while all the simulations overestimate the remaining indices (RTOT, R97.5pTOT, R20mm and Rx5day). The better performance of NEX than CORDEX in this regard may be because the NEX dataset is bias corrected while the CORDEX dataset is not. The bias-corrected CORDEX dataset was not available for the study.
In simulating the spatial distribution of the precipitation indices (Figure 4), NEX also performs better than CORDEX in most indices, except with regard to the extreme rainfall threshold (R97.5p). Consistent with Crétat et al. (2012), Figure 4 shows that the annual rainfall (RTOT) over the DMR increases from about 500 mm year −1 in the interior (in the west) to about 700 mm year −1 along the coast (in the east). While both datasets reproduce this east-west rainfall gradient and extend the maximum RTOT from the east to the northeastern top of the mountain range, their maximum RTOT is higher than observed. The bias in the NEX RTOT simulation (about 100 mm year −1 ) is lower than the uncertainty in the observation datasets (i.e., the standard deviation of the datasets; about 300 mm year −1 ), but that of CORDEX (more than 650 mm year −1 ) is more than the observation uncertainty. Also, the spatial correlation between NEX and the observation (r = 0.98) is higher than the one between CORDEX and the observation (r = 0.89). The CORDEX bias in RTOT can be linked to an error in the simulated number of WDAY ( Figure 4e) and in the amount of precipitation from extreme rainfall events (R97.5pTOT; Figure 4q). This suggests that the rainfall parameterization scheme in the CORDEX models may be too active in triggering or releasing precipitation over the mountain range. Several studies (Hohenegger et al., 2008;Ratna et al., 2014;Fosser et al., 2015;Favre et al., 2016;Pathak et al., 2019) have reported that many climate models overestimate rainfall over mountains because their convective parameterizations are too sensitive to the orographic lifting of moist air. Nevertheless, CORDEX does give an excellent representation of the extreme rainfall F I G U R E 4 Spatial distribution of the precipitation indices over the Drakensberg Mountain Range, as depicted by observation (Obs; the mean of ARC2-FEWS and CHIRPS) and simulation (i.e., the ensemble mean from Coordinated Regional Climate Downscaling Experiment (CORDEX) and National Aeronautics and Space Administration (NASA) Earth Exchange (NEX) Global Daily Downscaled Projections (NEX-GDDP) datasets. The difference and correlation between the two observation datasets are indicated with contours and correlation coefficients (r), respectively, in the Obs panels, while the model bias (with reference to Obs) and the correlation between simulation and observation are indicated in the other panels (CORDEX and NEX-GDDP). The precipitation indices are defined in Table 3 [Colour figure can be viewed at wileyonlinelibrary.com] threshold (R97.5p; r = .99) and performs better than NEX (r = .92) in this regard. The better performance of COR-DEX (than NEX) in R97.5p may be attributed to the physics in the CORDEX downscaling, which correctly captures the extreme rainfall events from the orographic lifting of moist air, whereas the NEX statistical downscaling does not incorporate such physics. This may also explain why NEX underestimates R97.5p over the DMR (by more than 8 mm day −1 on the mountain top). Nevertheless, while CORDEX outperforms NEX in simulating R97.5p, it underperforms NEX in simulating the other indices, because it produces too many extreme rainfall events and generates too much rainfall from these extreme events. Hertig et al. (2012) also showed that the frequency of extreme precipitation events over the Mediterranean mountains was better captured with statistical downscaling than with dynamical downscaling.
However, in general, our result corroborates previous studies that found no preference between RCM and Statistical Downscaling Model (SDM) in downscaling GCMs simulations over complex topography but advocated for a combination (hybrid) of both methods (Tang et al., 2016;Yhang et al., 2017;Kar et al., 2018;Tran Anh and Taniguchi, 2018). For instance, Yhang et al. (2017) found that while RCM overestimated precipitation over East Asia, SDM underestimated the precipitation variance, but a combination of both methods produced the best results in time and space. The hybrid dynamicalstatistical downscaling approach usually combines the ability of RCM to resolve fine-scale atmospheric features with the low computational cost of statistical downscaling (Tran Anh and Taniguchi, 2018).

| Projected changes in mean rainfall intensity
There are substantial differences between the CORDEX and NEX projections for RTOT over the DMR; the magnitude of the differences increases with the increase in the GWLs (Figures 5 and 6). For instance, with regard to the spatial average of the RTOT changes ( Figure 5), while the majority of the CORDEX simulations project a decrease in RTOT (the median is −2 mm year −1 at GWL1.5 and − 6 mm year −1 at GWL3.0), most NEX simulations suggest an increase (the median is 0.5 mm year −1 at GWL1.5 and 8 mm year −1 at GWL3.0), although the projections of both datasets are associated with large uncertainty. In terms of the spatial distribution of the RTOT changes (Figure 6), the greatest disagreement between the datasets occurs at high altitudes (>2000 m above sea level), where CORDEX projects a decrease in RTOT (up to 60 mm year −1 at GWL3.0), whereas NEX projects an increase (up to 30 mm year −1 at GWL3.0). At the lower altitudes (i.e., below 2000 m above sea level), both datasets project a decrease in RTOT, but the magnitude of the decrease is higher in the CORDEX projection than in the NEX. The discrepancy between CORDEX and NEX with regard to the RTOT projection may be attributed to a number of reasons. It may be due to the differences in the set of GCMs simulations downscaled by the two experiments. While CORDEX downscaled 11 GCMs and NEX downscaled 19 GCMs, and the two sets only have 9 GCMs in common. The discrepancy may also be due to how the two downscaling approaches handle climate feedback from local-scale forcing (land-sea boundary, topography, land use and land cover change). Several studies (e.g., Hewitson and Crane, 2006;Ayar et al., 2016) have shown that feedback from these local-scale features can modify (i.e., either amplify or reduce) large-scale climate signals (from GCMs) over an area. While the CORDEX RCMs used convection and turbulence parameterization schemes to capture the local-scale feedback, the NEX SDM does not have the schemes. Hence, the NEX SDM (which assumes a stationary relationship between GCM output and observation) may preserve the trends in the GCMs simulations, but the CORDEX RCMs may alter the trends because of the local-scale feedback.
Despite the disagreement between CORDEX and NEX on the RTOT projection, the two datasets agree on the WDAY and SDII projections (Figures 5 and 6). For instance, the ensemble means of both datasets project a decrease in WDAY and an increase in SDII over the DMR at all GWLs (suggesting fewer but more intense rainfall events in the future). This tendency for fewer but more intense rainfall events may be a regional extension to the Drakensberg area because previous studies (i.e., Abiodun et al., 2017;Pohl et al., 2017;Thoithi et al., 2021) have reported similar projections over most parts of southern Africa. Both datasets (NEX and CORDEX) agree that the magnitude of the projected changes in the indices grows as the GWL increases, except that, for WDAY, the growth is faster in CORDEX, and for SDII, the growth is faster in NEX. Furthermore, in both datasets, the WDAY and SDII projections are more robust than the RTOT projections as more than 70% of the simulations agree on the direction of the projected changes in WDAY and SDII.
Nevertheless, there are a few differences in the COR-DEX and NEX projections for these indices. For example, in the spatial distribution of the WDAY changes (Figure 6), F I G U R E 6 Projected changes in rainfall indices over southern Africa at four global warming levels (GWLs; GWL1.5, GWL2.0, GWL2.5 and GWL3.0), as depicted by Coordinated Regional Climate Downscaling Experiment (CORDEX) and National Aeronautics and Space Administration (NASA) Earth Exchange (NEX) datasets. The vertical strip (j) indicates where at least 80% of the simulations agree on the sign of the changes, while the horizontal strip (−) indicates where at least 80% of the simulations agree that the projected change is statistically significant (at 99% confidence level). The cross (+) shows where both conditions are satisfied; hence the change is robust [Colour figure can be viewed at wileyonlinelibrary.com] the magnitude of the projection is higher in CORDEX than in NEX (especially, in the area located more than 2,000 m above the sea level). In addition, going by the relationship between three variables (RTOT, WDAY and SDII), the consequence of RTOT and SDII changes on RTOT changes differ for both datasets. While CORDEX suggests that the annual water surplus from the increased SDII may be less than the annual water deficit from the reduced WDAYS (hence the decreased RTOT), NEX indicates that the annual water surplus from the increased SDII may be more than the annual water deficit from the decreased WDAYS, hence the increased RTOT. Nevertheless, Table 4 presents that projected changes in mean rainfall intensity over DMR are lower than the present-day inter-annual variability (i.e., standard deviation).

| Projected changes in duration indices
In both datasets (CORDEX and NEX), most simulations project a decrease in the mean wet spell duration (CWD) over the DMR and indicate that the magnitude of the projection varies non-linearly with the increase in GWLs ( Figure 5). However, the magnitude of the decrease is generally higher in the CORDEX dataset, where the projection uncertainty is also lower. The spatial distribution of projected changes is like that of RTOT ( Figure 6). While both datasets agree on a decreased CWD at the mountain base, they disagree on the changes near the mountain top, where CORDEX projects a decrease (>2 days) and NEX projects a weak increase (<0.5 day) or no change. This discrepancy, which is consistent with F I G U R E 6 (Continued) the one in RTOT, may be attributed to how the climate feedbacks from the topography are handled in the datasets. However, the wet spell metric used here needs to be considered with caution because it is very sensitive to threshold effects and is also less realistic under current climate conditions. In addition, the projected changes are lower than the present-day natural variability (Table 4).

| Projected changes in extreme rainfall indices
There is a good agreement between the CORDEX and NEX datasets with regard to the projected changes in the characteristics of extreme rainfall events over the DMR ( Figure 5). In both datasets, most simulations project an increase in mean extreme rainfall intensity (RA97.5p, RA97.5pTOT and Rx5days) and frequency (RA97.5pFREQ, R20mm and WEREFreq) over the DMR, showing that the magnitude of the increase generally grows with an increase in the GWLs ( Figure 5). For instance, with NEX, the median of the projected increase jumps from 1% (at GWL1.5) to 8% (at GWL3.0) in R95pTOT, from 2% (at GWL1.5) to 11% (at GWL3.0) in RA97.5pFREQ and from 2% (at GWL1.5) to 22% (at GWL3.0) in WEREFreq. Both datasets also agree that, at GWL3.0, the percentage increase in WEREFreq may be more than RA97.5pFREQ, suggesting that the future extreme rainfall event may be more from organized or mesoscale convective systems (MCS). However, given F I G U R E 6 (Continued) that the resolution of both datasets is too coarse to resolve MCS, more studies with higher resolution datasets are needed to investigate this further. The agreement between the datasets also features in the spatial distribution of the extreme precipitation indices (RA97.5p, RA97.5pTOT, RA97.5pFREQ and R20mm; Figure 6). For instance, in some indices (e.g., RA97.5p, RA97.5pTOT and RA97.5pFREQ), both datasets feature a region of increase that extends from the Indian Ocean to the mountain top, though the region of increase is wider in NEX than in CORDEX. The link between the increase of extreme rainfall activities over the Indian Ocean and the mountain top suggests a transport of warmer and moisture-laden (i.e., more buoyant) air from the Indian Ocean to the mountain top in the future. This is consistent with the findings from previous studies that linked precipitation over the DMR to moisture transport from the Indian Ocean (Tyson et al., 1976;Ndarana et al., 2021). For instance, Tyson et al. (1976) show that the north-eastern part of the gradient winds around the Indian Ocean high transports moist air to the Drakensberg.

| Scaling of extreme precipitation changes
There is an active debate in the climate community on how well the projected changes in extreme precipitation scale with Clausius-Clapeyron (CC) relationship (Allen F I G U R E 6 (Continued) and Ingram, 2002;Drobinski et al., 2018;Schröer and Kirchengast, 2018;Kendon et al., 2019). In light of the discussion, Figure 7 shows the scaling of extreme precipitation (R97.5p) changes with the regional temperature changes at various warming levels in comparison with the CC theoretical values (about 6.2% K −1 over DMR). The figure shows that, for both CORDEX and NEX, there is a large discrepancy among the simulations on the scaling. While some simulations suggest a higher scaling than the CC (i.e., super-CC scaling), some show a lower scaling (sub-CC scaling), and others feature a negative scaling (Figure 7a). However, the ensemble means of both datasets agree on a sub-CC scaling over DMR, with the maximum scaling (0.2 × CC) featuring over the eastern coast and the mountain's top. The projected sub-CC scaling of extreme precipitation changes over the DMR can be attributed to a number of factors. First, the projection assumes a constant relative humidity with the warming. Meanwhile, the changes in moisture availability may not increase as fast as temperature, leading to a decrease in relative humidity. This possibly explains why the magnitude of the scaling drops as the warming level increases (Figure 7). Using the dew point temperature (T d , which is the temperature to which an air parcel must be cooled to reach saturation) instead of temperature may give a more reliable scaling because T d is a measure of specific humidity translated to temperature using the CC relationship. Kendon et al. (2019) reported a super-CC scaling over most parts of Africa (especially over DMR) when T d was used for the scaling and a sub-CC scaling when temperature was used. Second, the projection also assumes that storm dynamics do not change with the warming. However, the environment in which the storms develop may change, altering characteristics (e.g., direction and vertical motion) of air that feeds the storm. Lastl, the convection parameterization schemes used in RCM and the empirical equations used in the SDM could also influence the scaling. Kendon et al. (2019) found higher scaling in convection-permitting simulations than in convection-parameterized simulations.

| The SOM classification of the future projections
To reveal the spatial patterns of the projection that might be hidden in the ensemble mean projections, this section discusses the SOM classification (4 × 3 nodes) of the CORDEX and NEX projections over the Drakensberg, focusing on the projected changes in the annual rainfall (RTOT; Figure 8), the extreme rainfall intensity (R97.5p; Figure 9) and the extreme rainfall frequency (R97.5pFREQ; Figure 10). The goal is to identify the most distinct or major patterns in the projected changes in these indices, and to examine how they differ from the ensemble mean patterns in Figure 6. Note that for any SOM classification, the most distinct patterns are usually the four edge nodes (i.e., Nodes 1, 4, 9 and 12).
For RTOT (Figure 8), the SOM analysis reveals that the spatial distribution of the projected changes can be grouped into three major patterns. The first pattern T A B L E 4 The simulated mean and SD (σ) of rainfall indices over DMR during the present-day and the projected changes under various global warming as depicted by the ensemble mean of CORDEX and NEX (in italics).  Table 2). Abbreviations: CORDEX, Coordinated Regional Climate Downscaling Experiment; DMR, Drakensberg Mountain Range; GCM, General circulation models; GWL, global warming levels; NEX, National Aeronautics and Space Administration (NASA) Earth Exchange.
features a decrease in RTOT (with different magnitudes) over the entire DMR domain (i.e., Figure 8a: Nodes 3,4,7,8,11 and 12). This most frequent pattern accounts for about 47% of the total RTOT projections (i.e., 50% in CORDEX and 43% in NEX). This implies that, while the ensemble mean of the NEX projections indicates an increase in RTOT over the domain, 43% of the projections actually indicate a decrease in RTOT over the domain. The second pattern shows an increase in RTOT (with different magnitudes) over most parts of the DMR domain (i.e., Figure 8: Nodes 5, 6, 9 and 10), accounting for 46% of the total RTOT projections (i.e., 36% in CORDEX and 57% in NEX). This also implies that, while the ensemble mean of the CORDEX projections indicates a decrease in RTOT over the domain (Figure 6), 36% of them suggest a decrease in RTOT over most parts of the domain. The third pattern is characterized by a band decrease in RTOT over the mountain top and along the coast, but an increase in RTOT elsewhere (Figure 8: Nodes 1 and 2).
This pattern, which only accounts for about 7.2% of the total RTOT projections (14% in CORDEX), does not feature in the NEX simulations. The SOM classification of the projected changes in R97.5p (Figure 9) can be generally grouped into two patterns. The first pattern features an increase in R97.5p over most parts of the DMR (i.e., Nodes 1, 2, 5, 11, 9 and 10) and accounts for 54% of the total R97.5p projections (50% in CORDEX and 58% in NEX). In contrast, the second pattern shows a decrease in R97.5p over most parts of the DMR (Nodes 3,4 7,8,11 and 12). While this pattern accounts for 46% of the total R97.5p projections (50% in CORDEX and 58% in NEX), it does not reflect in the ensemble mean projections for both NEX and CORDEX. The pattern of the SOM classification for R97.5pFREQ projections ( Figure 10) is similar to the pattern of R97.5p, and in most cases the simulations that projected an increase in R97.5p over most parts of the DMR also produced an increase in R97.5pFREQ over most parts of the F I G U R E 7 Scaling between future changes in extreme rainfall-temperature and the temperature over Drakensberg Mountain Range (DMR) at various Global Warming Levels as depicted by Coordinated Regional Climate Downscaling Experiment (CORDEX) and National Aeronautics and Space Administration (NASA) Earth Exchange (NEX) datasets. The scaling coefficient is obtained as the future changes in the logarithm of extreme precipitation intensity divided by future changes in mean temperature and has been normalized by the theoretical Clausius-Clapeyron (CC) relationship (about 6.0-6.2%K −1 over DMR); Hence, a value of 1 corresponds to CC scaling. The upper row shows the simulation spread of the scaling averaged over DMR, whereas the last two rows show the spatial distribution from the ensemble mean. The SOMs results present here can guide climate modelling communities in selecting a fewer but meaningful set of GCMs for downscaling (dynamical, statistical and hybrid) and for climate impact studies, with the aim of capturing the ensemble mean projection as well as the diversity in the projected changes using a subset of the GCMs. For example, Table 5 ranks the GCMs based on their contributions to the diverse patterns of changes in RTOT and R97.5p in the SOMs results.  has the topmost ranking because it features all the extreme patterns in the RTOT projections (i.e., Dipole, Wetter and Drier) as well as in the R97.5p projections (i.e., Stronger and Wetter) ( Table 5). In contrast, MRI-CGCM3 (w) has the lowest raking because it does not feature any extreme pattern in the RTOT or R97.5p projections. This implies that downscaling NorESM1-M (x) (instead of MRI-CGCM3 (w) ) projections would give more diverse patterns, hence wider range of uncertainty. Previous studies (e.g., Mizuta et al., 2014;Monerie et al., 2017) have employed similar approach (cluster analysis) to group GCMs based on model bias in historical climate simulations and on future climate projections (using temperature, precipitation and sea surface temperature). They also advocated for this type of ranking in a situation when it is not possible to downscale or use a large ensemble of models for impact studies.

| CONCLUSION
In this study, we investigated the projected changes in the characteristics of rainfall indices over the DMR at four different GWLs (1.5, 2.0, 2.5 and 3.0 C) under the RCP8.5 scenario. We compared statistical (NEX) and dynamical (CORDEX) downscaled projections of the CMIP5 GCMs. To investigate the capability of the two downscaled products to capture the observed characteristics of the rainfall indices over the DMR, we compared their biases in historical simulation of the past climate with observation uncertainties based on eight observation datasets (ARC, CHIRPS, PERSIANN, TAMSAT, AgCFSR, AgMERRA, WFDEI-CRU and WFDEI-GPCC). The model ensemble spatial means were compared using box plots, whereas the projected change in the spatial distribution of rainfall indices was compared using SOMs. The results of the study can be summarized as follows: • The eight observation datasets agree on the characteristics of precipitation (intensity-frequency curves and precipitation indices) over the DMR, albeit with some discrepancies. • The CORDEX and NEX datasets give a realistic simulation of precipitation characteristics over the DMR and, in most cases, their biases lie within the observation uncertainty. • NEX performs better than CORDEX in reproducing most of the precipitation characteristics (e.g., intensityfrequency curves, RTOT, WDAY, SDII), but CORDEX performs better than NEX in simulating the threshold of extreme rainfall (R97.5p). • The ensemble means of both datasets agree on a future increase of SDII, RA97.5p, RA97.5pTOT, RA97.5p-FREQ and WEREFreq and on a decrease of WDAYS and CWD, but they disagree on the projected changes of RTOT, for which CORDEX projects an increase over most parts of the domain, while NEX indicates a decrease. • The SOM analysis, which reveals projected patterns that are invisible in the ensemble mean results, shows the four most important projections with a combination of annual rainfall and extreme precipitation events (intensity and frequency): (a) increase in both annual precipitation and extreme precipitation events; (b) decrease in both annual precipitation and extreme precipitation events; (c) increase in annual precipitation but decrease in extreme precipitation events; and (d) decrease in annual precipitation but increase in extreme precipitation events.
The results of this study can provide a basis for developing climate change adaptation and mitigation strategies over the Drakensberg. Based on the results of the SOMs, Table 6 suggests the most possible three scenarios regarding the changes in mean precipitation intensity and extreme rainfall events over the DMR domain, indicates the potential impacts, and suggests adaptation strategies. For example, the projected increase in extreme precipitation events over the DMR suggests an increased future risk of floods and geotechnical hazards. Hunter et al. (2016) reported that an extreme rainfall event over the Main Ranges of the Canadian Rocky Mountains and over the Bow River watershed (Canada) in June 2013 produced flooding that led to geotechnical and hydrologic hazards. The occurrence of such geotechnical hazards over the DMR may be detrimental to the lives and socioeconomic activities of the communities at the foot of the DMR. The projected decrease in annual precipitation may furthermore lead to a future increase in droughts in the DMR. Droughts in this area could have severe implications for water resources and disaster management. It may result in a drop in river flow level in rivers such as T A B L E 6 Potential impacts of the projected changes in mean rainfall intensity and extreme rainfall events over the DMR and the proposed adaptation strategies.

Projection
Potential impacts Adaptation strategies Increase in both annual rainfall and extreme rainfall events • Increased risk of floods, landslides, geotechnical and hydrologic hazards • Degradation of ecological services, which may involve risks to lives, livelihoods and lifestyles the Tugela and Vaal, which rise from the upper reaches of the Drakensberg. These rivers provide water for agriculture, hydro-electric power generation and home use downstream. In addition, drought may negatively affect the Lesotho Highlands Water Project, which is a water supply project with a hydropower component, developed in partnership between the governments of Lesotho and South Africa. Putting an adequate plan or policy in place to mitigate these impacts would help to reduce the associated risks.
However, to provide more robust information that is need for policy makers, the results of this study can be improved in several ways. For example, the projected droughts over the Drakensberg from the SOM analysis are based on precipitation alone; they do not however account for the role of potential evapotranspiration. Meanwhile, Abiodun et al. (2019) showed that accounting for atmospheric evaporative demand (potential evapotranspiration) in the drought projections (i.e., using the Standardized Precipitation Evapotranspiration Index [SPEI]) produced more severe droughts than using only precipitation (or the Standardized Precipitation Index) for such projections. Hence, future studies could consider incorporating the influence of potential evapotranspiration on future drought projections over this mountain range. In addition, the downscaled datasets do not have all the dynamical variables needed to investigate the dynamics behind the projected changes. Future studies could investigate the dynamics behind the projected changes using a series of sensitivity experiments with regional climate models (forced by GCM simulations) or a global model with a variable horizontal grid resolution. Furthermore, the present study has used regional climate simulations in which convection is activated through physical parameterizations of convective processes. Several studies have shown that explicit convection may produce sensibly different results and changes in rainfall properties (e.g., Kendon et al., 2019). Hence, analysing simulations with high resolution and explicit convection for the projection would provide more detailed information over a complex topography like the Drakensberg. However, given that it is practically impossible to use such a high-resolution model to downscale all the GCM simulations analysed in the present study, the ranking of the GCM simulations (Table 5) could guide in the optimal selection of the GCM simulations for the downscaling.