Identiﬁcation and simulation of space–time variability of past hydrological drought events in the Limpopo River basin, southern Africa

. Droughts are widespread natural hazards and in many regions their frequency seems to be increasing. A ﬁner-resolution version (0.05 ◦ × 0.05 ◦ ) of the continental-scale hydrological model PCRaster Global Water Balance (PCR-GLOBWB) was set up for the Limpopo River basin, one of the most water-stressed basins on the African continent. An irrigation module was included to account for large irrigated areas of the basin. The ﬁner resolution model was used to analyse hydrological droughts in the Limpopo River basin in the period 1979–2010 with a view to identifying severe droughts that have occurred in the basin. Evaporation, soil moisture, groundwater storage and runoff estimates from the model were derived at a spatial resolution of 0.05 ◦ (ap-proximately 5 km) on a daily timescale for the entire basin. PCR-GLOBWB was forced with daily precipitation and temperature obtained from the ERA-Interim global atmospheric reanalysis product from the European Centre for Medium-Range Weather Forecasts. Two agricultural drought indicators were computed: the Evapotranspiration Deﬁcit Index (ETDI) and the Root Stress Anomaly Index (RSAI). Hydrological drought was characterised using the Standardized Runoff Index (SRI) and the Groundwater Resource Index (GRI), which make use of the streamﬂow and groundwater storage resulting from the model. Other more widely used meteorological drought indicators, such as the Standardized Precipitation Index (SPI) and the Standardized Precipitation Evaporation Index (SPEI), were also computed for different aggregation periods. Results show that a carefully set-up, process-based model that makes use of the best available input data can identify hydrological droughts even if the model is largely uncalibrated. The indicators considered are able to represent the most severe droughts in the basin and to some extent identify the spatial variability of droughts. More-over, results show the importance of computing indicators that can be related to hydrological droughts, and how these add value to the identiﬁcation of hydrological droughts and ﬂoods and the temporal evolution of events that would otherwise not have been apparent when considering only meteorological indicators. In some cases, meteorological indicators alone fail to capture the severity of the hydrological drought. Therefore, a combination of some of these indicators (e.g. SPEI-3, SRI-6 and SPI-12 computed together) is found to be a useful measure for identifying agricultural to long-term hydrological droughts in the Limpopo River basin. Additionally, it was possible to undertake a characterisation of the drought severity in the basin, indicated by its time of occurrence, duration and intensity.


P. Trambauer et al.: Identification and simulation of past hydrological droughts in the Limpopo River basin
result of climate variability and climate change (IPCC, 2007;Patz et al., 2005;Sheffield and Wood, 2008;Lehner et al., 2006). Moreover, and probably more importantly, the rapid increase in world population will certainly aggravate water shortage on local and regional scales. The study of droughts and drought management planning has received increasing attention in recent years as a consequence.
Drought monitoring is a key step in drought management, requiring appropriate indicators to be defined by which different types of drought can be identified. Meteorological, agricultural and hydrological drought indicators are available to characterise different types of droughts. The best known indicators are the Standardized Precipitation Index (SPI; McKee et al., 1993) and the Palmer Drought Severity Index (PDSI, Palmer, 1965;Dube and Sekhwela, 2007;Alley, 1984); both are primarily meteorological drought indices. The SPI uses only precipitation in its computation, and the PDSI uses precipitation, soil moisture and temperature. However, the timescale of drought that the PDSI addresses is often not clear (Keyantash and Dracup, 2002) and will usually be determined by the timescale of the data set; Vicente-Serrano et al. (2010b) indicate that the monthly PDSI is generally correlated with the Standardized Precipitation Evaporation Index (SPEI) at timescales of about 9-12 months. While the computation of the PDSI is complex, applied to a fixed time window and difficult to interpret, the SPI is easy to compute and to interpret in a probabilistic sense, is spatially invariant and can be tailored to a time window appropriate to a user's interest (Guttman, 1998). Alley (1984) and Vicente-Serrano et al. (2010b) also highlight several limitations of the PDSI, such as not allowing for the distinction of different types of drought (i.e. hydrological, meteorological and agricultural) as it has a fixed temporal scale. The PDSI has other derivatives such as the Palmer Hydrological Drought Index (PHDI) for hydrological long-term droughts, Palmer Z Index for short-term monthly agricultural droughts and the Crop Moisture Index (CMI) for short-term weekly agricultural droughts. The empirical PDSI method developed in the United States, is still widely used in the US but is gradually being substituted by other indicators in other regions (Keyantash and Dracup, 2002) as a result of its limitations. The SPI can be computed for different timescales by accumulating the precipitation time series over the time period of interest (typically 3 months for the SPI-3, 6 months for the SPI-6, and 12 months for the . The SPI has shown to be highly correlated with indicators of agricultural drought, hydrological drought and groundwater drought. The SPI-3 has a high temporal variability that is associated with shortto-medium range meteorological anomalies that can result in anomalous soil moisture and crop evolution, and it can therefore be used as an indication of agricultural drought. The SPI-6 has a higher correlation with hydrological droughts, mainly represented by low anomalies in runoff. The SPI-12 and SPI-24 have a lower temporal variability and point to major and long-duration drought events whose impacts may extend to groundwater. The widely used SPI does, however, have its limitations, mainly because it is based only on precipitation data. An extension of the SPI was proposed by Vicente-Serrano et al. (2010b), called the Standardized Precipitation Evaporation Index (SPEI), which is based on precipitation and potential evaporation. In a way, it combines the sensitivity of the PDSI to changes in evaporation demand with the capacity of the SPI to represent droughts on multitemporal scales (Vicente-Serrano et al., 2010b).
Together with the development of the first drought indicators, hydrological models were also used for agricultural and hydrological drought assessment. Schulze (1984) applied the Agricultural Catchments Research Unit (ACRU) hydrological model in Natal, South Africa, to compare the severity of the 1979-1983 drought with other drought events in the previous 50 years. He identified hydrological modelling as a potentially powerful tool in drought assessment. Moreover, he indicated that it is necessary to distinguish between different types of droughts, as droughts in terms of water resources do not necessarily coincide with droughts from the crop production point of view.
In recent years, several new indicators have been developed to characterise the different types of drought. Although drought indicators are mostly used to characterise past droughts and monitor current droughts, forecasting of these indicators at different spatial and temporal scales is gaining considerable attention.
In this study we extend a continental-scale framework for drought forecasting in Africa, which is currently under development (Barbosa et al., 2013), and apply this to the Limpopo Basin in southern Africa, one of the most waterstressed basins in Africa. The Limpopo River basin is expected to face even more serious water scarcity issues in the future, limiting economic development in the basin (Zhu and Ringler, 2012). To apply this framework at the regional scale, a finer-resolution version of the global hydrological model PCRaster Global Water Balance (PCR-GLOBWB) was adapted to regional conditions in the basin. We model hydrological droughts and their space-time variability using a process-based distributed hydrological model in the (semi-) arid Limpopo Basin. The model was tested by comparing the simulated hydrological and agricultural drought indicators in the period 1979-2010 with reported historic drought events in the same period. We derive a number of different drought indicators from the model results (see Table 1), such as the ETDI (Evapotranspiration Deficit Index; Narasimhan and Srinivasan, 2005), the RSAI (Root Stress Anomaly Index), the SRI (Standardized Runoff Index; Shukla and Wood, 2008) and the GRI (Groundwater Resource Index; Mendicino et al., 2008). While the SRI is based on river discharge at a particular river section, the ETDI, RSAI and GRI are spatial indicators that can be estimated for any location in the basin. The ETDI and RSAI are directly related to water availability for vegetation with or without irrigation, and the GRI is related to groundwater storage. Moreover, we compute the widely known meteorological drought indicators SPI and SPEI at different aggregation periods to verify the correlation of the different aggregation periods for these indices and the different types of droughts. Table 1 presents the derived indicators with a description of the purpose and the type of drought each indicator represents. The aim of this study is to assess the ability of different drought indicators to reconstruct the history of droughts in a highly water-stressed, semi-arid basin. Moreover, we investigate whether widely used meteorological indicators for drought identification can be complemented with indicators that incorporate hydrological processes.

Study area: Limpopo River basin
The Limpopo River basin has a drainage area of approximately 415 000 km 2 and is shared by four countries: South Africa (45 %), Botswana (20 %), Mozambique (20 %) and Zimbabwe (15 %) (Fig. 1). The climate in the basin ranges from tropical dry savannah and hot dry steppe to cool temperatures in the mountainous regions. Although a large part of the basin is located in a semi-arid area, the upper part of the basin is located in the Kalahari Desert, where it is particularly arid. Aridity, however, decreases further downstream. Rainfall in the basin is characterised as being seasonal and unreliable, causing frequent droughts, but floods can also occur in the rainy season. The average annual rainfall in the basin is approximately 530 mm year −1 , ranging from 200 to 1200 mm year −1 , and occurs mainly in the summer months (October to April) (LBPTC, 2010). Arid and semi-arid regions are generally characterised by low and erratic rainfall, high interannual rainfall variability and a low rainfall-to-potential-evaporation ratio. This leads to the ratio of runoff to rainfall being low on the annual scale. Hydrological modelling possesses considerable challenges in such a region. A detailed discussion on problems related to rainfall-runoff modelling in arid and semi-arid regions can be found in Pilgrim et al. (1988).
The runoff coefficient (RC = runoff / precipitation) of the Limpopo Basin is remarkably low. For the station at Chokwe (no. 24), which is the station with the largest drainage area among the discharge stations available in this study (Fig. 1), the runoff coefficient is just 4.3 % for the naturalised discharge and a mere 1.7 % for the observed discharge (without naturalisation). Note that the naturalised discharge is estimated as observed discharge plus the estimated abstractions. These runoff coefficients are strikingly low: out of 539 mm year −1 of annual rainfall only 23 mm year −1 (basin average) turns into runoff annually, including abstraction. This means that even a small error in estimates of precipitation and evaporation could result in a large error in the runoff. Moreover, the uncertainty in the rainfall input could easily be larger than the runoff coefficient (4.3 %) of the basin. Runoff coefficients for other selected stations in the basin (highlighted in Fig. 1, right panel) are presented in Table 2.
The basin is also highly modified, as is evident from the observed and naturalised runoff. This adds an additional challenge to modelling this basin. For example, for the largest drainage outlet available (no. 24), the observed annual discharge is only some 39 % of the naturalised discharge, which means that the abstractions in the basin amount to 61 % of the total runoff. Irrigation water demand takes up the largest share. The total estimated present demand in the basin is about 4700 × 10 6 m 3 year −1 . The total natural runoff generated from rainfall is approximately 7,200 × 10 6 m 3 year −1 , showing that a significant portion of the runoff generated in the basin is currently used.

Data for the hydrological model
The digital elevation model (DEM) we used is based on the Hydro1k Africa (USGS EROS, 2006). The majority of the parameters (maps) required for the model (soil layer depths, soil storage capacity, hydraulic conductivity, etc.) were derived mainly from three maps and their derived properties: the Digital Soil Map of the World (FAO, 2003), the distribution of vegetation types from Global Land Cover Characterization (GLCC) (USGS EROS, 2002;Hagemann, 2002) and the lithological map of the world (Dürr et al., 2005). From the soil map, 73 different soil types were distinguished in the basin. The irrigated area was obtained from the global map of irrigated areas in 5 arcmin resolution based on Siebert et al. (2007) and FAO (1997). We computed the monthly irrigation intensities per grid cell using the irrigated area map, the irrigation water requirement data per riparian country in the basin and the irrigation cropping pattern zones (FAO, 1997). All meteorological forcing data used (precipitation, daily temperature, daily minimum and maximum temperature at 2 m) are the same as in Trambauer et al. (2014) and are based on the ERA-Interim (ERAI, Dee et al., 2011) reanalysis data set from the European Centre for Medium-Range Weather Forecasts (ECMWF). This data set covers the period from January 1979 to the present day with a horizontal resolution of approximately 0.7 • and 62 vertical levels. A comprehensive description of the ERAI product is available in Dee et al. (2011). The ERA-Interim precipitation data used with the present model were corrected with GPCP v2.1 (product of the Global Precipitation Climatology Project) to reduce the bias with measured products (Balsamo et al., 2010). The GPCP v2.1 data are the monthly climatology provided globally at a 2.5 • × 2.5 • resolution, covering the period from 1979 to September 2009. The data set combines the precipitation information available from several sources (satellite data, rain gauge data, etc.) into a merged product (Huffman et al., 2009;Szczypta et al., 2011). From September 2009 to December 2010, the mean monthly ERAI precipitation was corrected using a mean bias coefficient based on the climatology of the bias correction coefficients used for the period 1979-2009. While this only corrects for systematic biases, this was the only option available at the time, as a new version of GPCP (version 2.2) was not available. Temperature data is used for the computation of the reference potential evaporation needed to force the hydrological model. In this study the Hargreaves formula was used. This method uses only temperature data (minimum, maximum and average), so it requires less parameterisation than Penman-Monteith, with the disadvantage that it is less sensitive to climatic input data, with a possibly reduction of dynamics and accuracy. However, it leads to a notably smaller sensitivity to error in climatic inputs (Hargreaves and Allen, 2003). Moreover, the potential evaporation derived from the Penman-Monteith equation and Hargreaves equation result in very similar values throughout Africa, and the choice of the method used for the computation of the reference potential evaporation appears to have minor effects on the results of the actual evaporation for southern Africa (Trambauer et al., 2014). For this study, the ERAI data were obtained for the period of 1979-2010. These were converted to the same spatial resolution as the continental-scale model using bilinear interpolation to downscale from the ERAI grid to the 0.5 • model grid. ERAI is archived using an irregular grid (reduced Gaussian) over the domain and thus an interpolation was inevitable to be able to use it in the model.
Runoff data were obtained from the Global Runoff Data Centre (GRDC; http://grdc.bafg.de/), the Department of Water Affairs in the Republic of South Africa and ARA-Sul (Administração Regional de Aguas do Sul, Mozambique). Runoff stations that had data available up until recent years, with relatively few missing data, are presented in Fig. 1. Most of these stations are in the South African part of the basin as almost no data could be found from stations in the other riparian countries. The subbasins draining to each hydrometric station are named after the station number.

General description
A process-based distributed hydrological (water balance) model based on PCR-GLOBWB ( van Beek and Bierkens, 2009) is used. First the global-scale model was adapted to the continent of Africa (Trambauer et al., 2014). A higherresolution version (0.05 • × 0.05 • ) of the continental model (0.5 • × 0.5 • ) was applied for the Limpopo River basin. The PCR-GLOBWB was one of the 16 different land surface and hydrological models reviewed (Trambauer et al., 2013), and it was identified as one of the hydrological models that can potentially be used for hydrological drought studies in large river basins in Africa. PCR-GLOBWB is in many ways similar to other global hydrological models, but it has many improved features, such as improved schemes for sub-grid parameterisation of surface runoff, interflow and baseflow, a kinematic-wave-based routing for the surface water flow, dynamic inundation of floodplains, and a reservoir scheme (van Beek and Bierkens, 2009; van Beek, 2008). On a cell-by-cell basis and at a daily time step, the model computes the water storage in two vertically stacked soil layers (max. depth 0.3 and 1.2 m) and an underlying groundwater layer, as well as computing the water exchange between the layers and between the top layer and the atmosphere. It also calculates canopy interception and snow storage. Within a grid cell, the sub-grid variability is taken into account considering tall and short vegetation, open water and different soil types. Crop factors are specified on a monthly basis for short-and tall-vegetation fractions, as well as for the open-water fraction within each cell. These crop factors are calculated as a function of the leaf area index (LAI) as well as of the crop factors for bare soil and under fullcover conditions (van Beek et al., 2011). Monthly climatology of LAI is estimated for each GLCC (Global Land Cover Characterization)-type, using LAI values per type for dormancy and growing season from Hagemann et al. (1999). LAI is then used to compute the crop factor per vegetation type according to the FAO guidelines (Allen et al., 1998). The total specific runoff of a cell consists of the surface runoff (saturation excess), snowmelt runoff (after infiltration), interflow (from the second soil layer) and baseflow (from the lowest reservoir as groundwater). River discharge is calculated by accumulating and routing specific runoff along the drainage network and including dynamic storage effects and evaporative losses from lakes and wetlands (van Beek and Bierkens, 2009;van Beek, 2008). The default PCR-GLOBWB model does not explicitly consider irrigated areas but the version of the model used here includes an irrigation module to account for the highly modified hydrological processes in the irrigated areas of the basin.

Drought indicators
The meteorological drought indicators used in this study are computed only from meteorological variables: precipitation and potential evaporation. Agricultural and hydrological indicators, on the other hand, are computed from the results of the hydrological model and therefore account for effects of soil, land use, groundwater characteristics, etc., in the basin. The indicators used in this study are described below.

Meteorological drought indicators Standardized Precipitation Index (SPI)
The SPI was developed by McKee et al. (1993) and it interprets rainfall as a standardised departure with respect to a rainfall probability distribution. It requires fitting the precipitation time series to a gamma distribution function, which is then transformed to a normal distribution allowing the comparison between different locations. The SPI [−] is then computed as the discrete precipitation anomaly of the transformed data divided by the standard deviation of the transformed data (Keyantash and Dracup, 2002;McKee et al., 1993). SPI values mainly range from −2 (extremely dry) to 2 (extremely wet).

Standardized Precipitation Evaporation Index (SPEI)
Instead of using only precipitation as in the SPI, the SPEI uses the difference between precipitation (P ) and potential evaporation (PET), i.e. D = P − PET, and the PET is computed following the Thornthwaite method (Vicente-Serrano et al., 2010a, b). The calculated D values are aggregated at different timescales, following the same procedure as for the SPI. A log-logistic probability function is then fitted to the data series of D, and the function is then standardised following the classical approximation of Abramowitz and Stegun (1965). The SPEI also ranges between −2 and 2; the average value of the SPEI is 0, and the standard deviation is 1 (Vicente-Serrano et al., 2010a, b).

Agricultural drought indicators
Agricultural droughts are defined as the lack of soil moisture to fulfil crop demands, and therefore the agriculture sector is normally the first to be affected by a drought. In this study we characterise agricultural droughts by means of two spatially distributed indicators defined as described in the following.

Evapotranspiration Deficit Index (ETDI)
The ETDI (Narasimhan and Srinivasan, 2005) is computed from the anomaly of water stress to its long-term average. The monthly water stress ratio (WS [0-1]) is computed as: where PET and AET are the monthly reference potential evaporation and monthly actual evaporation, respectively. The monthly water stress anomaly (WSA) is calculated as where MWS y,m is the long-term median of water stress of month m, maxMWS m is the long-term maximum water stress of month m, minWS m is the long-term minimum water stress of month m, and WS y,m is the monthly water stress ratio (y = 1979-2010 and m = 1-12). Narasimhan and Srinivasan (2005) scaled the ETDI to between −4 and 4 to be comparable with the PDSI. Here, we used the same scaling procedure but amended this to scale the ETDI to between −2 and 2 to make it comparable to the SPI, SPEI and SRI: ETDI y,m = 0.5 ETDI y,m−1 + WSA y,m 100 . (4)

Root Stress Anomaly Index (RSAI)
The "root stress" (RS) is a spatial indicator of the available soil moisture, or the lack of it, in the root zone. The root stress varies from 0 to 1, where 0 indicates that the soil water availability in the root zone is at field capacity and 1 indicates that the soil water availability in the root zone is zero and the plant is under maximum water stress. The RSAI is computed similarly to the ETDI described above. The monthly root stress anomaly (RSA) is calculated as where MRS m is the long-term median root stress of month m, maxMRS m is the long-term maximum root stress of month m, minRS m is the long-term minimum root stress of month m, and RS y,m is the monthly root stress (y = 1979-2010 and m = 1-12). The root stress anomaly index, scaled to between −2 and 2 (using the same procedure as Narasimhan and Srinivasan, 2005) is

Hydrological drought indicators
For the characterisation of hydrological droughts we used the commonly applied Standardized Runoff Index (SRI; Shukla and Wood, 2008) for streamflow and the Groundwater Resource Index (GRI; Mendicino et al., 2008) for groundwater storage.

Standardized Runoff Index (SRI)
The SRI follows the same concept as the SPI and is defined as a "unit standard normal deviate associated with the percentile Hydrol. Earth Syst. Sci., 18, 2925-2942, 2014 www.hydrol-earth-syst-sci.net/18/2925/2014/ of hydrologic runoff accumulated over a specific duration" (Shukla and Wood, 2008). To compute the SRI the simulated runoff time series is fitted to a probability density function (a gamma distribution is used here), and the function is used to estimate the cumulative probability of the runoff of interest for a specific month and temporal scale. The cumulative probability is then transformed to the standardised normal distribution with a mean of 0 and a variance of 1 (Shukla and Wood, 2008).

Groundwater Resource Index (GRI)
The GRI y,m is suggested as a standardisation of the monthly values of groundwater storage (detention) without any transformation (Mendicino et al., 2008): where S y,m is the value of the groundwater storage for the year y and the month m, and µ S,m and σ S,m are respectively the mean and the standard deviation of the groundwater storage S simulated for the month m in a defined number of years (32 years in this case). The same classification that is used for the SPI (between −2 and 2) is applied to the GRI (Wanders et al., 2010).

Identification of past droughts and primary characterisation of drought severity
To identify past droughts, the drought indicators described were calculated for the period 1979-2010 for the Limpopo River basin, resulting in times series of monthly indicator maps. The maps allow for the visualisation of the spatial variability of the indicators in the basin. The SPI, SPEI and SRI were computed for different aggregation periods (1, 3, 6, 12 and 24). All the indicators were then aggregated over several subbasins resulting in times series for each indicator. The historical subbasin-averaged indicators were then compared. Maps of the indicators are also compared for specific years to show the spatial variability of the indicators and the extent of the droughts. All indices considered were scaled to range between −2 and 2. Based on the SPI values, droughts may be classified into mild (0 > SPI ≥ −1), moderate (−1 > SPI ≥ −1.5), severe (−1.5 > SPI ≥ −2) and extreme (SPI < −2) (Lloyd-Hughes and Saunders, 2002; see Table 3). For the SPI and SPEI, the spatially averaged indicators are no longer related to a probability of occurrence. However, we still use the same thresholds for the characterisation of the subbasin aggregated droughts, as we understand that the resulting indicators would not be very different from the computation of these indicators with aggregated precipitation and potential evaporation. For agricultural (ETDI and RSAI) and groundwater indicators (GRI) this is not the case as these are not defined based on a probability of occurrence.
where Iv is the indicator value; j starts with the first month of a drought and continues to increase until the end of the drought (x) for any of the i timescales. The DS [months] would be numerically equivalent to the drought duration if the drought had an intensity (value) of −1.0 for each month (McKee et al., 1993).

Hydrological model performance
Given the complexity of the basin for hydrological modelling, particularly due to the arid or semi-arid nature, the model results are quite satisfactory, especially for the larger subbasins. Runoff estimates from the hydrological model were verified with observed runoff on a monthly basis. For a number of the runoff stations tested, the coefficient of determination (R 2 ) values varied from about 0.45 to as good as 0.92. In a review of model application and evaluation, Moriasi et al. (2007) recommended three quantitative statistics for model evaluation: Nash-Sutcliffe efficiency (NSE), percent bias (PBIAS) and the ratio of the root mean square error to the standard deviation of the measured data (RSR). They also specified ranges for these statistics for a "satisfactory" model performance (NSE > 0.5, RSR ≤ 0.70 and PBIAS ± 25 % for streamflow). However, PBIAS is highly influenced by uncertainty in the observed data (Moriasi et al., 2007). Given the potential problems in observed flow data in South Africa, reported by the Water Research Commission (2009), such as poor accuracy of the rating table, particularly at low flows, and the inability to measure high flows, we do not evaluate our results based on PBIAS. The evaluation measures NSE and RSR together with the coefficient of  Table 4. We do not calibrate parameters based on these evaluation measures, but we use them as a simple test of concordance. Based on the ranges proposed by Moriasi et al. (2007), the model performance is found to be satisfactory for four out of six runoff stations.  sons 1982sons /1983sons , 1986sons /1987sons , 1991sons /1992sons and 1994sons /1995sons . The drought of 1991sons /1992 was the most severe in the region in recent history. After the year 2000, important droughts include the years 2002/2003/2004 and 2005/2006. Droughts in the Limpopo River basin also show significant spatial variability. A study covering only the Botswana part of the basin documents a severe drought that occurred in 1984 (Dube and Sekhwela, 2007). However, in that year no documentation of drought in the other parts of the basin was found. Figure 2 presents the RSAI and ETDI for the most severe drought in recent history (1991/1992), for the very dry year 1982/1983, for a wet year (1999/2000) and for a year with both dry and wet conditions at different locations in the basin (1984/1985). The geographic variability of the RSAI seems to be slightly higher than that of the ETDI. These indicators provide information for the assessment of agricultural droughts. The figure shows that both indicators, computed from different outputs of the hydrological model (actual evaporation and soil moisture), produce similar results and are able to reproduce the dry or wet conditions in the basin. This is also supported by Fig. 3, which shows the fraction of the Limpopo Basin under moderate to extreme agricultural drought, i.e. Iv < −1.0. Both indicators illustrate that a large part of the basin was under at least moderate agricultural drought conditions for the years with recorded drought events. Figure 4 shows the SRI values (1, 3, 6 and 12 months) from 1980 to 2010 computed from the simulated runoff at station 24. The dotted grey line at the threshold value of −1 is used to identify moderate droughts, with the moderate drought considered to start when the indicator falls below the threshold, and stop when the indicator goes above the threshold. The simulated SRI clearly identifies the severe hydrological droughts of 1982/1983 and 1991/1992   with missing data and the computation of the SRI requires a monthly runoff data set for a continuous period without missing data. The Groundwater Resource Index (GRI) presented in Fig. 5 for the same selected years shows the years 1991/1992 and 1982/1983 to be drier than normal, but the intensity of the drought appears to be quite low (not severe). The year 1984/1985, selected as it presents both dry and wet conditions at different locations in the basin, does not show this spatial variability for the GRI. This was to be expected, due to the persistence of the groundwater storage and low intensity of indicators of drought or wetness in this year in different locations of the basin. The intensity of the extremely wet year 1999/2000 is well represented, suggesting that the GRI is skewed. This is likely due to the fact that the GRI is not transformed into the normal space. Moreover, the distribution of values is constrained by the capacity of the groundwater reservoir in the hydrological model. Mendicino et al. (2008) applied this indicator in a Mediterranean climate but the skewness test of normality showed that their series from January to September were normally distributed, while the series of October to December were not normally distributed. However, they indicate that the values of groundwater storage in the last winter months and in spring were more important. For this indicator to be applied independently of the climate and basin conditions, it should probably be transformed into the normal space.

Comparison of drought indicators
The computed indicators were averaged for the whole basin as well as for the selected subbasins. Time series of the resulting indicators were compared for the whole 1980-2010 period. Figure 6 presents the time series of aggregated drought indicators for subbasin 24. Note that the subbasins are named after the hydrometric station number. Figure 6 compares the agricultural, hydrological and groundwater drought indicators. The agricultural indicators ETDI and RSAI are compared with the meteorological drought indicators SPI and SPEI with the short aggregation period (3 months) that is commonly used as indicators of agricultural droughts. Figure 6 (upper plot) shows that the indices are mostly in phase, correctly representing the occurrence of dry and wet years, and the intensities of the events are in general quite similar. The hydrological drought indicator SRI-6 is compared with the meteorological drought indicators SPI-6 and SPEI-6 (upper middle plot). All three indicators follow roughly the same pattern, but the fluctuation of the SRI seems to be slightly lower than that of the meteorological indices (SPI and SPEI). This is probably due to the higher persistence of streamflow when compared to precipitation. Moreover, it is clearly visible from Fig. 6 that the temporal variability or fluctuation of the indicators reduces when moving from drought indicators associated with agricultural drought to those associated with hydrological drought. This means that several mild agricultural droughts do not progress further to hydrological droughts. Moreover, to identify groundwater droughts, or major drought events, the time series of the GRI is compared to the time series of meteorological and hydrological drought indicators with long aggregation periods (SPI-12, SPEI-12, SRI-12, SPI-24, SPEI-24, SRI-24) (see Fig. 6, lower middle and lower plots). The plots show that as the variability of the indicator reduces further, the number of multi-year, prolonged droughts increases. However, for groundwater droughts, only two events (1982/1983 and 1991/1992) are identified as moderate to severe droughts (Iv < −1). The plots again show that, in general, the temporal variability of the runoff-derived indicator (SRI) is lower than that of the meteorological indicators (SPI and SPEI). The GRI shows much less temporal variability than the other indices and does not identify any extreme events, with the exception of the flood of 1999/2000. Similar results using the GRI were found by Wanders et al. (2010), who indicate that the GRI has a very low number of droughts with a high average duration. Moreover, a study of Peters and Van Lanen (2003) investigated groundwater droughts for two climatically contrasting regimes. For the semi-arid regime they found multi-annual droughts to occur frequently. They indicate that the effect of the groundwater system is to pool erratically occurring dry months into prolonged groundwater droughts for the semi-arid climate. Table 5 presents a correlation matrix between all the indicators considered in this study for subbasin 24. Similar correlation results were found for the other subbasins. The table shows that the agricultural drought indicators ETDI and RSAI have the highest correlation with the SPEI-3, SPEI-6, SPI-3, SPI-6 and with the SRI with low aggregation periods (1 to 3 months). For every station the correlation between the agricultural indicators and the SPEI is slightly higher than with the SPI. While the hydrological drought indicators SRI-6 and SRI-12 present the highest correlation with the meteorological drought indicators SPI-12 and SPEI-12, the extended hydrological drought indicator SRI-24 is better correlated with the meteorological drought indicators SPI-24 and SPEI-24. The GRI shows the highest correlation with the SRI-6 and SRI-12. This makes sense, given the direct connection between groundwater and runoff, where groundwater (baseflow) contributes to the total runoff. Figures 7-9 present the monthly spatial mean time series of drought indicators for subbasins 1, 18 and 20, respectively. The averaged indicators for subbasins 24 and 1, the two largest subbasins considered, are almost identical (see Figs. 6 and 7). Figure 8 shows that even though the general pattern of the time series for subbasin 18 is similar to that found for subbasins 24 and 1, some differences are noticeable. For example, Fig. 8 shows a clear drought period for subbasin 18 in the years 1984/1985/1986, which is not apparent for subbasins 24 and 1. These localised drought events that affected the upper part of the basin were not apparent for the lower part of the basin. This was also observed in Fig. 2. Moreover, the extreme floods that occurred in the lower part of the basin in 1999/2000 are much less severe in the upstream parts of the basin. For example, Fig. 9 shows that for subbasin 20 (the smallest subbasin considered), the flood of 1996/1997 was more severe than that of 1999/2000. Similarly, while the drought of 2003/2004 is quite mild when averaged over the largest selected subbasin (no. 24), it is quite severe for subbasin 20 (similar to the droughts of 1983/1984 and 1991/1992).
For the four subbasins a short but intense agricultural drought is noticeable at the beginning of the 2005/2006 season, but this did not progress to an extended hydrological drought. This is consistent with the literature, which indicates that this season was delayed and after a dry start to the season, good rainfall occurred from the second half of December (Department of Agriculture of South Africa, 2006). In subbasin 18 (Fig. 8), even though meteorological indicators (SPI-6, SPEI-6, SPI-12 and SPEI-12) suggest that the 1986/1987 season was near normal to wet, the hydrological indicators (SRI-6, SRI-12) point to a dry runoff year. Measured runoff at this station indicates that the year 1986/1987 was indeed a dry year. This seems similar to what was found by Peters and van Lanen (2003); for longer aggregations periods an accumulation of successive short anomalies can lead to an overall hydrological drought. Similarly, meteorological indicators suggest that the floods of 1996/1997 and 1999/2000 in the lower part of the basin were of a similar magnitude. However, records indicate that the flood of 1999/2000 was much more extreme than the one of 1996/1997(WMO, 2012. This can be seen clearly in the hydrological drought indicators SRI-6, SRI-12 and SRI-24. The help identify the spatial and temporal evolution of drought and flood events that would otherwise not have been apparent when considering only meteorological indicators. We also computed drought severities (DS (months)) resulting from the different indicators as explained in Sect. 3.3 (Eq. 9). The droughts of 1982/1983, 1986/1987, 1991/1992, 1994/1995, 2002/2003/2004 and 2005/2006 are identified as being among the most severe droughts, but the end month of these drought events varies for the different indicators. The indicators with higher aggregation periods (e.g. 12 and 24 months), which have a lower temporal variability, generally point to longer droughts (multi-year droughts) with higher persistence than indicators with lower aggregation periods (agricultural droughts). For example, while the agricultural indicators suggest that the extreme drought of 1991/1992 was over by the end of 1992 or beginning of 1993, the indicators that represent hydrological droughts signal that this drought only ended at the end of 1993. Moreover, for the SRI-12, GRI, SPI-24 and SPEI-24, this multiyear drought lasts until 1994/1995. As an example for subbasin 24, Fig. 10 presents the duration and severity of the six most severe recorded droughts as identified by the meteorological drought indicator SPEI aggregated for different periods to represent agricultural, hydrological and extended hydrological droughts (multi-year droughts). The graph shows that the multi-year droughts resulting from the accumulation of shorter successive droughts are the most severe as a result of the duration. These droughts can be the most hazardous, as a succession of mild droughts that can initially seem nonproblematic can result in very severe droughts if they last for P. Trambauer et al.: Identification and simulation of past hydrological droughts in the Limpopo River basin a long time. The average intensity of these droughts is generally lower than that of the agricultural droughts, which can be very intense but often of shorter duration.

Conclusions
Very low runoff coefficients and high rainfall variability pose major challenges in modelling hydrological droughts in (semi-)arid basins. Small errors in the meteorological forcing and estimation of evaporation may result in significant errors in the runoff estimation. This also implies that model calibration, if any, should be applied cautiously to maintain the physical meaning of model parameters. We opted to apply a process-based model and parameterise it on the basis of the best available input data without additional calibration.
In the process we ensured that we were using reliable data sets and interpolated or aggregated them with care to prepare spatially distributed parameter maps. The model is able to simulate hydrological drought-related indices reasonably well. We have derived a number of different drought indicators from the model results, such as the ETDI, RSAI, SRI and GRI. While the SRI is based on river runoff at a particular river section, the ETDI, RSAI and GRI are spatial indicators that can be estimated at any location in the basin. The ETDI and RSAI are directly related to water availability for vegetation with or without irrigation, and the GRI is related to groundwater storage. Moreover, we computed the widely known drought indicators SPI and SPEI at different aggregation periods to verify the correlation of the  different aggregation periods for these indicators and the different types of droughts. All the indicators considered (with the exception of the GRI) are able to represent the most severe droughts in the basin and to identify the spatial variability of the droughts. Our results show that even though meteorological indicators with different aggregation periods serve to characterise droughts reasonably well, there is added value in computing indicators based on the hydrological model for the identification of hydrological droughts or floods and their severity. The indicators also help identify the spatial and temporal evolution of drought and flood events that would otherwise not have been apparent when considering only meteorological indicators.
The RSAI follows the ETDI to a great extent, and the ETDI is quite well represented by the SPEI-3 and the SPEI-6. This indicates that in the absence of actual evaporation and soil moisture data which are required to compute the ETDI and RSAI, the meteorological indicator SPEI-3, which considers both precipitation and potential evaporation and is reasonably easy to compute, may be used as an indicator of agricultural droughts. For discharge we observe some added value in computing the SRI. Even though the SPI can give a reasonable indication of drought conditions, computing the SRI can be more effective for the identification of hydrological droughts. The groundwater indicator GRI mostly remains near normal conditions. A combination of different indicators, such as the SPEI-3, SRI-6 and SPI-12 (computed together), can be an effective way to characterise agricultural to long-term hydrological droughts in the Limpopo River basin.