The Benefit of Continuous Hydrological Modelling for Drought Hazard Assessment in Small and Coastal Ungauged Basins: A Case Study in Southern Italy

Abstract

Rainfall-runoff modelling in small and ungauged basins represents one of the most common practices in hydrology. However, it remains a challenging task for researchers and practitioners, in particular in a climate change context and in areas subject to drought risk. When discharge observations are not available, empirical or event-based approaches are commonly used. However, these schemes can be affected by several relevant assumptions. In the last years, continuous models have been developed in order to address the major drawbacks of event-based approaches. With this goal in mind, in this work we applied a synthetic rainfall generation model (STORAGE; stochastic rainfall generator), constituting the implementation of a modified version of Neymann-Scott rectangular pulse (NSRP) model, and a continuous rainfall-runoff framework (COSMO4SUB; continuous simulation modelling for small and ungauged basins) specifically designed for ungauged basins within a climate change context. The modeling approach allows one to investigate the drought hazard using specific indicators for rainfall and runoff in a small watershed located in southern Italy. Results show that the investigated area seems to tend to a mild/moderate drought in a future time period of approximately 30 years, with a decrease in seasonal water volumes availability in the range of 15–30%. Finally, our results confirm that the continuous modelling is suitable for rapid and effective design simulations supporting drought hazard assessment.

1. Introduction

Global climate change has had a growing impact on water resources management and drought hazard assessment, and the world population has paid the cost of meteorological hazards [1,2].

In 2015, the United Nations encouraged economic, environmental, and social advancements in order to obtain a sustainable future for the entire world population (in particular in a climate change context), introducing the Sustainable Development Goals (SDGs) [3]. Integrated water resources management is part of the SDGs, but today only 2.5% of all water on Earth can be used for human domestic needs [4]. This aspect as increased researcher and practitioners’ interest on water scarcity and droughts.

From several studies, it appears evident that is not possible to provide a unique drought definition since this is related to hydro-meteorological conditions, geographical locations, and water demands in different regions around the world [5]. Hence, the drought definition can be classified using four categories [6,7]: meteorological [8,9,10,11,12], hydrological [13,14], agricultural [15], and socio-economic [16]. The first three categories provide direct methods to quantify drought as a physical phenomenon, while the last one measures its direct effect on the social production and on the human life.

The complexity of the processes related to drought makes difficult the identification of a unique diagnostic criterion, which should be able to define the status and trend of the phenomenon. Among the various methods developed by the scientific community, the application of suitable “indicators” seems to have an undoubted usefulness. Particularly interesting is the streamflow drought index (SDI) method, used to predict the drought onset and its duration. This method uses the cumulative flow rate of the river as well as the spatial extent of drought for assessing and monitoring the hydrological droughts [17]. When it is possible to integrate rainfall and vegetation data, it seems more appropriate to use the normalised difference vegetation index (NDVI) and standard precipitation index (SPI) [18]. Also, the Palmer drought severity index (PDSI) has been largely used when temperature and precipitation data are available. However, monthly PDSI values do not capture droughts on time scales lower than 12 months [19].

Focusing in the Mediterranean area, it is not easy to select, develop, and apply indicators to assess the drought hazard, since the region is particularly vulnerable due to climatic and geomorphological conditions and to consolidated anthropic pressure. In this context, desertification research projects such as MEDALUS can be particularly important, aiming to implement ‘Plans National Action’, referred to the ESAs model (environmental sensitive areas to desertification) [20,21]. MEDALUS provided a methodology to standardize all possible indicators in order to identify and discriminate regions at high risk of desertification in the Mediterranean area. In detail, in the MEDALUS model, the vulnerability is substantially connected to climatic, morphodynamic, pedological, vegetational, and anthropic factors. This method was widely applied to several case studies in northern Africa, in the Middle East, and in other Mediterranean-type ecosystems [22]. In Italy, a long-time interval study of vulnerable land growth has been conducted using the ESA approach [23]. The ESA output, called ESAI (environmental sensitivity area index), is an indicator system which examines more than 10 variables assessing climate, soil, vegetation, and intensity in land use [24].

Essentially, the main problem in defining indicators is the need for a robust time series [25,26,27]. In particular, the World Meteorological Organization (WMO) suggested to consider at least 30 years of continuous rainfall data related to the resolution of interest for specific applications [26]. Such a circumstance is not easy to obtain in small and ungauged basins [28,29], making the estimation of runoff consequent to rainfall difficult [30,31,32]. In this context, the use of stochastic rainfall generators (SRGs) could help in the analysis of rainfall processes [33]. A SRG generally presents a simple mathematical formulation and low computational costs, and it can be easily used for simulating long rainfall time series by assuming any selected climate change scenario (or eventually in a stationary context).

Moreover, the availability of a long rainfall time series could allow one to apply a continuous rainfall-runoff transformation, and consequently to obtain the whole runoff time series for the selected basin. In doing so, we could have available a design simulation, allowing the robust estimation of drought characteristics and indexes.

The aim of the manuscript is hence to demonstrate the importance of continuous hydrological modelling for the drought hazard assessment in ungauged basins. In order to do that, in this work we first applied two recently developed hydrological models. The former is a synthetic rainfall generation model (STORAGE) [34,35] that was applied in a climate change scenario. The latter is a continuous rainfall-runoff model (COSMO4SUB) [36,37,38]. After that, we examined several hydrological drought indices based on the modelled rainfall and runoff time series in order to provide an in-depth analysis for a small basin located in southern Italy (Apulia region) that is potentially subject to drought. This approach, never used in the selected study area, represents an important novelty. In fact, the Apulia region unfortunately has a scarce availability of runoff data, making any type of hydrological study complicated. The present manuscript is organized as follows. In Section 2, the materials and methods are described. In particular, the selected case study, the proposed SRG, the continuous rainfall-runoff model, and the selected indicators are descripted. In Section 3, our results are presented and discussed. In Section 4, the conclusions are reported.

2. Materials and Methods

2.1. Case Study Description

The study area pertains to the hydrographic basin of the Cillarese river (Figure 1), a torrential watercourse which has historically represented a fundamental element for the socio-economic development of the area. In particular, the area is included within the municipalities of Oria, Torre Santa Susanna, Latiano, Mesagne and Brindisi, and belongs to the province of Brindisi, in the Salento peninsula, a large area in the southern part of the Apulia region (southern Italy). This region is historically prone both to flooding and to drought [39,40,41], and one of the known major problems of the area is the overexploitation of the groundwater [42,43,44]. The analyzed river originates in the municipal territory of Torre Santa Susanna at 103 m above sea level, and extends for about 28 km until it reaches the harbor of Brindisi city in the Adriatic Sea.

In 1980, a dam was built and the course of the Cillarese river in the bottom valley was blocked, allowing for the formation of a water reservoir of over 4 million m³ (with a surface area equal to 2.76 km²) for the water supply of the Brindisi industrial development area.

Selecting the dam location as the outlet river cross section, the Cillarese hydrographic basin is characterized by an extension of about 132 km² and an average slope of 1.3%.

Following the construction of the reservoir on the Cillarese river, a fauna protection oasis (extension equal to 1.27 km²) was implemented [45]. Furthermore, the area is classified as SIR (Site of Regional Interest n. IT9140012), in compliance with the European Commission Habitats Directive—92/43/EEC [46]. Based on the ESAI, the area is classified at maximum risk of desertification (category C3—ESAI > 1.53) and, therefore, is among the protected areas of the Apulia Region [47]. These data are confirmed by the Apulia Region River Basin Authority that declared that, on average, the whole Apulia region accounts for only 500 mm of rainfall with 60–80 wet days per year (in Italy, a day is defined as wet if the daily rainfall depth is equal or greater than 1 mm) [48].

Concerning the available rainfall data, the rain gauge located in Brindisi (10 km far from the catchment’s centroid) was considered as representative for the Cillarese basin. The time series of annual maximum rainfall (AMR) for durations equal to 1, 3, 6, 12, and 24 h (observation period: 1936–2019), and daily data (observation period: 1950–2019) are available and were used for calibration of the adopted synthetic rainfall generator. Data were provided by Apulia Region Civil Protection.

For the selected case study, the rainfall data present a mean annual precipitation (MAP) of 612 mm and a mean annual number of wet days equal to 69. Concerning mean values of seasonal precipitation, we have 203 mm for December-January-February (DJF), 134 mm for March-April-May (MAM), 60 mm for June-July-August (JJA), and 215 mm for September-October-November (SON).

2.2. The Synthetic Rainfall Generation Model

In this study the synthetic rainfall generator implemented in STORAGE software [34,35] has been used. In detail, the model constitutes a modified version of the Neyman–Scott Rectangular Pulses (NSRP) model [49,50]. In the basic version of the NSRP model, five quantities, which are considered as random variables (and hence follow assigned probability distributions) play a crucial role:

• The waiting time (assumed as exponentially distributed) between the occurrences of two consecutive storms;
• The number of rain cells (also named as bursts or pulses) in each storm. This quantity is considered as a geometric random variable;
• The waiting time (assumed as exponentially distributed) between the occurrences of a storm and of an associated cell;
• The intensity and the duration (both considered as exponentially distributed) of each cell inside a storm.

STORAGE is characterized by two main innovations as respect to the usual NSRP model described in literature. First, the parametric estimation is carried out by using data series (i.e., annual maxima rainfall, annual and monthly cumulative rainfall, annual number of wet days) which are usually longer than observed high-resolution series (which are commonly used for the calibration of a usual NSRP model). Second, the seasonality is modelled using series of goniometric (i.e., sine and cosine) functions. This approach makes STORAGE less demanding, in terms of the total needed number of parameters, as it usually occurs for basic versions of NSRP models.

In this work, a transient version of STORAGE was adopted [34], in which some parameter trends are assumed (see Section 3) on the basis of the available projections of regional climate models (RCMs) for the investigated area, in order to obtain a high-resolution synthetic rainfall time series which respects the trends hypothesized by RCMs at coarser scales. The hypothesis that the statistics of a really-averaged rainfall simulated by climate models reflect changes at the point rain gauge scale is commonly used in literature [34].

In detail, the publication of the Italian Institute for Environmental Protection and Research (ISPRA) [51] considered four RCMs (ALADIN52, CMCC-CCLM4, GUF-CCLM4, and LMDZ4), and for each one we have two scenarios of representative concentration pathways (RCP) of total radiative forcing (i.e., the cumulative measure of human emissions of greenhouse gasses from all sources expressed in W/m²): RCP 4.5 W/m² (intermediate emissions) and RCP 8.5 W/m² (high emissions). With respect to the reference period 1971–2000, the results of the RCM simulations for the spatial cells comprising Cillarese basin and related to the future period from 2000 to 2070 are:

• A decrease of the annual cumulative precipitation value, which is comprised between −100 and −50 mm;
• A modest increase (no greater than 5 mm) for the annual maximum daily rainfall.

In detail, the application of STORAGE here performed consisted of the following steps. First, under the hypothesis of stationary process, STORAGE was calibrated according to [35], i.e., by employing the parameters of amount-duration-frequency (ADF) curves, mean values for MAP, mean annual number of wet days, and mean values of DJF, MAM, JJA, and SON rainfall depths.

Second, the stationary version of STORAGE was validated by generating 500-year time series at 5-min resolution and by comparing the synthetic and observed data in terms of frequency distributions for AMR, annual and seasonal rainfall, and annual number of wet days.

Third, the climate change analysis was carried out by assuming some trends for STORAGE parameters, using those that provided compatible results with RCM scenarios in terms of variation for maximum daily rainfall and annual precipitation. In detail, the adopted scenario combines the following hypotheses:

• A linear increasing trend of about 26% in 50 years concerning the mean waiting time between two consecutive storms;
• A linear increasing trend of about 22% in 50 years concerning the mean value of intensity of rainfall cells;
• A linear decreasing trend of about 14% in 50 years concerning the mean value of duration of rainfall cells.

Fourth, 500 synthetic time series, each one concerning 51 years of continuous 5-min rainfall depths, were generated. The first generated year, denoted as “0”, has the features of the calibrated stationary process, while the successive years are characterized by the assumed parameter trends. In this case, analysis with multiple realizations is necessary as the process is not assumed as stationary in this step, and then investigation of only one long time series is not possible, since the ergodicity property cannot be considered [52]. The generated 500 rainfall realizations were the input data for the COSMO4SUB input data, and the results have been averaged on the 500 rainfall realizations as explained in the following paragraph.

2.3. The Continuous Rainfall-Runoff Model

The continuous rainfall-runoff model used in the present manuscript is the COSMO4SUB framework introduced by Grimaldi et al. [36,37] and recently updated by Grimaldi et al. [38]. COSMO4SUB is characterized by four steps: (1) the rainfall scenario definition; (2) the excess rainfall estimation; (3) the excess rainfall-runoff transformation; and (4) the design simulation strategy.

Regarding step (1), the synthetic rainfall time series (generated by STORAGE) were used as input of COSMO4SUB.

Regarding step (2), COSMO4SUB applies the Curve Number (CN) for Green-Ampt (GA)-CN4GA-procedure [53]. CN4GA is a mixed model that first applies the National Resources Conservation Service (NRCS) CN method [54] in order to transform, at the event scale, the gross rainfall cumulative depth in excess rainfall depth (i.e., determining the cumulative infiltration depth), and then uses the GA equation [55] to distribute in time, within the rainfall event, the excess rainfall. The GA equation parameters are automatically estimated by CN4GA constraining the GA equation to furnish the same total excess rainfall depth and the same ponding time equal to that estimated by the CN method. It is noteworthy that, in order to adapt at the continuous modelling approach, the CN4GA procedure that was originally developed at event scale (i.e., a sort of rainfall event identification named separation time, Ts) is needed. Basically, the separation time represents the dry (or almost dry) period that should be spent waiting so that the soil initial abstraction became effective again. In detail, model parameters inherent to this step are essentially CN (which can vary from 0, i.e., all rainfall becomes infiltration and no runoff is formed, to 100 (i.e., all rainfall becomes runoff and no infiltration is allowed)), which can be estimated based on land cover and soil type, and Ts, which must be assigned by user. In the present work, following the suggestions of Grimaldi et al. [38], we considered 24 h as value for Ts, while CN was estimated as 70.4 in AMC-II condition (i.e., in soil moisture average condition). Starting from the CN value in AMC-II condition, and based on the cumulative value of rainfall occurred in the 5 antecedent days period, the CN value is automatically adjusted by COSMO4SUB for the generic rainfall event assigning a higher or lower value in case of AMC-III (wet soil) or AMC-I (dry soil) condition, according to the original NRCS-CN formulation [54].

Regarding step (3), the width function based instantaneous unit hydrograph (WFIUH) model based on the partial contributing area concept [56] is applied for convoluting excess rainfall into runoff time series. Regarding the WFIUH, this work has used the approach proposed by Petroselli [56], dividing the basin area in two parts. The first part takes into account the basin saturated area, responsible of the basin “fast” response to rainfall (the surface runoff). The second part takes into account the basin unsaturated area, responsible of the basin “slow” response to rainfall (the subsurface runoff). The division of the basin in the two saturated and unsaturated areas (mutually exclusive) is performed following the steps introduced by Petroselli [56] using the topographic wetness index and assuming a percentage of basin saturated area as a function of the cumulative rainfall value (characterizing the single rainfall event in the generated rainfall time series). In doing so, such WFIUH optimizes the available DEM, land cover and soil type information, and it estimates the travel time distribution of watershed DEM cells thanks to: (i) determining automatically the flow paths; (ii) distinguishing between hillslope-channel; and (iii) determining the distributed hillslope surface (in the saturated area) and subsurface (in the unsaturated area) flow velocities based on local slope, land cover, and soil data. Concerning the channel cells, the river network velocity is automatically calibrated assuming that the projection on the time axis of the WFIUH (i.e., its median point) is equal to the basin lag time, which is assumed as the 60% of the concentration time [54,56]. In this step, the only model parameter is basically Tc, determining the river network velocity, while local slope, land cover, and soil data determine the hillslope cells velocities.

Regarding step (4), and based on what previously expressed, the 500 rainfall realizations generated by STORAGE were fed as input data in COSMO4SUB, allowing one to obtain the corresponding 500 runoff scenarios (each one constituting of 51 years of runoff at 5-min resolution) and to identify the metrics described in the following paragraph.

2.4. The Investigated Metrics Based on the Generated Rainfall and Runoff Time Series

Based on the 500 generated rainfall realizations and the consequent 500 runoff time series, the following metrics have been investigated. These indicators have been chosen since they are particularly efficient in individuating environmentally sensitive areas prone to desertification [20,21].

Regarding the rainfall time series, we make reference to some of the metrics present in the European Climate Assessment & Dataset (https://www.ecad.eu//indicesextremes/indicesdictionary.php (accessed on 31 January 2022)), and in particular:

• MAP: yearly cumulative precipitation amount (mm);
• WD1: yearly cumulative number of wet days (amount of precipitation greater or equal than 1 mm) (days);
• DRYD: yearly cumulative number of dry days (amount of precipitation lower than 1 mm) (days);
• GSTP: yearly cumulative growing season (from April to October) precipitation (mm);
• NGSTP: yearly cumulative non-growing season precipitation (from November to March) (mm);
• DP10: yearly cumulative number of days where the daily precipitation amount is greater or equal than 10 mm (days);

The previous rainfall metrics have been averaged on the 500 rainfall realizations, each year.

Regarding the runoff time series, we analyze:

• TNGR: total number of isolated (Ts = 24 h) rainfall events in the 51 years for each generated rainfall time series (-);
• TNER: total number of isolated excess rainfall events in the 51 years for each generated rainfall time series (-);
• CUMVOL: yearly cumulated volume, averaged on the 500 rainfall realizations (m³);
• SDI: streamflow drought index (-), i.e., a well-known index using monthly streamflow values (here averaged on the 500 rainfall realizations) and a process of normalization associated for developing a drought index based upon streamflow data [48]. Literature states that for SDI < −2.0 there is an extreme drought condition, for −2.0 < SDI < −1.5 there is a severe drought condition, for −1.5 < SDI < −1.0 there is a moderate drought condition, for −1.0 < SDI < 0 there is a mild drought condition, and for SDI > 0 there is a non-drought condition;
• FDC: flow duration curves (m³/month). Flow duration curves (here averaged on the 500 rainfall realizations) represent cumulative frequency curves that show the amount of time when specified volumes are equaled or exceeded during a given period.

The analysis of the previously described metrics allows performing an in-depth analysis for the selected case study in terms of drought hazard.

3. Results and Discussion

3.1. Analysis of the Synthetic Rainfall Time Series

As stated, STORAGE was first calibrated and validated under the hypothesis of stationary process. The results of this first step of rainfall analysis are reported on EV1 (Extreme Value distribution of type 1, also named as Gumbel distribution) probabilistic plots (Figure 2) for AMR series and on Gaussian plots (Figure 3) for the other investigated series. The fitting of the sample distributions can be considered as acceptable, except for 3-h AMR (where we have a marked over-estimation), annual number of wet days (marked under-estimation) and MAM (marked under-estimation).

Then, the transient analysis was carried out by assuming the previously discussed parameters trends. In particular, for each metric it was possible to estimate a trend, as defined in Figure 4, on the time series of 51 years.

The results are illustrated in Figure 4, where six plots show the variability of the six investigated indices. Regarding DRYD, the trend shows that the index is growing, confirming that the study area has a propensity to drought, according also to what indicated by Apulia Region River Basin Authority [48]. For the same reason, it is observed a reduction of the values related to the metrics WD1 and DP10. The obtained results are in line with previous literature studies, which show that in many not-drought regions decreases of dry days and increases of wet days have been observed [57,58].

Regarding the rainfall volumetric metrics, it is interesting to note that the three considered indices (CUMVOL, SDI, FDC) present a decreasing trend. However, the GSTP index, evaluated as yearly cumulative growing season (from April to October) precipitation, shows an almost constant trend, with small fluctuations compared to the average value of the considered time sample. Zhao et al. [59] reported similar GSTP values in a high agricultural drought risk area. In general, also comparing with ecological study based on GSTP [60], our results show an inclination to drought. Analysing the scientific literature, MAP is an indicator largely used to define plant species that could grow and survive in arid and semiarid areas that, in a climate change context, could become drier and warmer. In particular, we can define for MAP a threshold between “low” and “high” MAP environments. Some authors define low-MAP environments the areas characterized by a MAP value lower than 400 mm [61].

In the present study, MAP is always greater than 500 mm, so the investigated area could not be classified as having high risk of drought. The study of Stuart-Haëntjens et al. [62] predicted forest and grasslands resistance and resilience to extreme drought using a MAP threshold equal to 828 mm. So, we can state that in the investigated area the presence of forests and grasslands would be at risk. However, MAP results showed that the annual precipitation decreased during the observed period, confirming some of the results pointed out by Buttafuoco et al. [63] in Southern Italy, and in our opinion confirming that the selected study area is prone to drought risk.

3.2. Results of the Continuous Rainfall-Runoff Modeling

Due to the assumed scenario of rainfall change during the years, we show in Figure 5 the total number of isolated rainfall events (Ts = 24 h) and the corresponding total number of excess rainfall events occurring in the 500 generated rainfall time series. This preliminary analysis shows that the 500 generated rainfall time series present a stable number of rainfall events, which ranges between 2125 and 2425, with an average number of 2304 (i.e., approximately 45 rainfall events each year). Regarding the total number of excess rainfall events, this is mainly due to the study area CN value, equal to 70.4 in AMC-II condition, but also to the shifts in AMC-I (dry soil) or AMC-III (wet soil) condition that can occur based on the rainfall temporal distribution of the generic rainfall time series realization.

As can be seen from Figure 5, the total number of excess rainfall events occurring in the 500 generated rainfall series also present a stable value ranging between 217 and 333, with an average number of 278 (i.e., approximately 5.4 excess rainfall events each year). Such value could seem a low value, but it can be attributed to the not so high CN value, which allows the majority of rainfall to infiltrate and hence to produce groundwater recharge and baseflow as respect to surface runoff. Indeed, as reported in the regional water protection master plan [47], the study area is characterized by high permeability of the karstic substrate that favors infiltration of rainfall, and there are few surface rivers concentrated in the northern area (Candelaro, Cervaro, Carapelle, and Ofanto). The scarcity of superficial water courses has generated in the past in the Puglia’s aquifers important overexploitation (especially by agriculture) [43].

In any case, in our opinion the limited variability within the 500 generated rainfall time series of the total number of rainfall events and excess rainfall events justifies the following analysis of the runoff time series.

Regarding the runoff analysis, in Figure 6 we show the seasonal (DJF, MAM, JJA, SON) time series of the cumulated yearly volumes, averaged on the 500 rainfall realizations. As can be seen from Figure 6, and as expected due to the Italian hydrological regime, the greater number of volumes is concentrated in the DJF and SON period, while the lower amount is related to the summer period, being approximately the 40% of the volume pertaining to the winter period. Except for the summer period, where the cumulative runoff appears to be quite constant, in the other three seasons we can observe a decrease in the water availability. In particular, in the 50-year time period investigated, the cumulative runoff decreases are approximately 30% in the DJF season, 33% in the MAM season, and 15% in the SON season.

In Figure 7 we show the yearly SDI time series, averaged on the 500 rainfall realizations and calculated based on runoff monthly values. As it can be seen from Figure 7, the adopted rainfall scenario tends to generate, in average, a mild drought at the end of the investigated 50 years, but fluctuations in the trend, producing a moderate drought, are visible after approximately 35 years.

Finally, in Figure 8, we concentrate on the FDC, represented in terms of cumulative runoff that is equaled or exceeded for a certain number of months each year. FDC have been calculated for each year in descending order for the cumulative monthly values of runoff. Then, for each month the 500 realizations are averaged. For the sake of simplicity, in Figure 8, we show the average of the FDC related to the first 10 years and the last 10 years of the investigated period. The decreasing trend in water availability is evident looking from Figure 8, with decreases of approximately 20% for the first two months and 12% for the remaining 10 months.

In general, the obtained results are in line with the trends available in the technical literature in the same area. For example, during the 1951–2000 period, total annual precipitation across Puglia region has significantly decreased at a rate of 23.9 mm/decade, which could lead to a strong reduction of the mean value if this trend would continue for one century more. Generally, long term observations from meteorological stations in the Puglia region show trends towards warmer and drier conditions during the second half of the 20th Century [64]. The trend analysis of considered metrics indicates an increasing drought in the case study area, confirming previous studies [64,65].

4. Conclusions

In this paper, a continuous rainfall-runoff modeling is proposed to identify the hydrological response in an ungauged coastal area prone to drought hazard. The methodology consists of two recently developed simple and parsimonious models, namely STORAGE and COSMO4SUB. The first is a synthetic rainfall generation model and the second is a continuous rainfall-runoff model. Projected time series of rainfall and runoff have been analyzed in a climate change context, selecting appropriate metrics for drought analysis. Results show that the investigated area seems to tend to a mild/moderate drought in a time period of approximately 30 years, and that a decrease in water volumes availability in the range 15–30% can be expected. In conclusion, in our opinion, our results confirm that the continuous modelling approach presented here is suitable for a rapid and effective design simulation supporting drought hazard assessment. The use of STORAGE and COSMO4SUB seems to represent an easy modeling approach that could be employed for evaluating possible measures needed to minimize the impact of drought (for instance, the reuse of agricultural wastewaters that is proposed by several researchers [66]) but that today is not so diffused in the Apulia Region [43,67].