Drought in the Twenty ‐ First Century in a Water ‐ Rich Region: Modeling Study of the Wabash River Watershed, USA

: Climate change is expected to alter drought regimes across North America throughout the twenty ‐ first century, and, consequently, future drought risk may not resemble the past. To explore the implications of nonstationary drought risk, this study combined a calibrated, regional ‐ scale hydrological model with statistically downscaled climate projections and standardized drought indices to identify intra ‐ annual patterns in the response of meteorological, soil moisture, and hydrological drought to climate change. We focus on a historically water ‐ rich, highly agricultural watershed in the US Midwest—the Wabash River Basin. The results show likely increases in the frequency of soil moisture and hydrological drought, despite minimal changes in the frequency of meteorological drought. We use multiple linear regression models to interpret these results in the context of climate warming and show that increasing temperatures amplify soil moisture and hydrological drought, with the same amount of precipitation yielding significantly lower soil moisture and significantly lower runoff in the future than in the past. The novel methodology presented in this study can be transferred to other regions and used to understand how the relationship between meteorological drought and soil moisture/hydrological drought will change under continued climate warming.


Introduction
Droughts-defined as prolonged periods of below-average water availability [1]-are one of the main types of weather-related disasters. Drought can be classified into three main typesmeteorological drought, soil moisture drought, and hydrological drought-representing belownormal conditions in precipitation, soil moisture, and streamflow, respectively. Meteorological droughts often precede other types of drought; however, the common impacts of drought, including reduced agricultural yield, forest fires, and water scarcity, are directly related to soil moisture and hydrological drought and only indirectly related to meteorological drought [2,3]. For example, the 1988 and 2012 droughts in the continental United States (US) resulted in an estimated US $40 billion and US $30 billion in mostly agricultural losses [4].
Due to the growing world population and increasing water demands, the adverse impacts of droughts are likely to worsen in the future. Recent drought studies [5][6][7] have shown an increasing trend in drought extent and affected population, and climate change is predicted to lead to more extremes [8]. Due to these factors (increasing population + increasing water demand + climate change), research on the climate change impacts on drought are economically and socially necessary.
Many previous studies have investigated climate change impacts on drought in North America [7,9,10] and regional studies have focused on the western, central, and eastern United States, see [11][12][13] among others. Most of these studies focus on meteorological and soil moisture drought, and a common conclusion for water-rich regions, like the US Midwest, is that soil moisture droughts will increase in frequency and severity despite decreasing frequency of meteorological droughts [13,14]. The contrasting changes in meteorological drought and soil moisture drought highlight the importance of incorporating temperature into assessments of future drought regimes. While the role of temperature in changing drought regimes has received increasing attention over the last decade [15,16], the implications of a nonstationary climate on future drought risk are not yet well understood [17].
To explore the potential implications of nonstationary drought risk, this study focuses on a historically water-rich, highly agricultural region-the Wabash River Basin. The Wabash Basin covers more than 65% of the state of Indiana, and more than 72% of the watershed is agricultural, with 66% cultivated crops and 6% pasture/hay. Cropland in this region is primarily rainfed [18] and thus exhibits relatively high drought sensitivity [19]. This study uses a large ensemble of 61 hydrological model simulations to investigate climate change impacts on seasonal water availability in the Wabash River Basin. Specifically, this study focuses on seasonal water availability and addresses three main questions: (1) How will climate change impact the annual and intra-annual water balance? (2) Will there be significant changes in the frequency and severity of short-term droughts and do the changes vary seasonally? and (3) Will the relationships between meteorological drought and soil moisture and hydrological drought change with continued climate warming?
The remainder of this paper is organized as follows. Section 2 presents the region of study, the model setup, calibration, and validation, along with methods used to calculate drought indices. The results are presented in Section 3 and discussed in Section 4. Conclusions are provided in Section 5.

Study Area
The domain for this study is the Wabash River Basin located in Midwest USA, which is the largest northern tributary of the Ohio River. The Wabash Basin drains over 90,000 km 2 , including most of the state of Indiana, part of eastern Illinois, and a small section of western Ohio ( Figure 1). The watershed has a humid continental climate, with a mean annual temperature of 11.3 °C and mean annual precipitation of 122.7 cm over the 1971-2000 period.

Hydrologic Model
The Soil and Water Assessment Tool (SWAT) hydrological model (version 2012) was used to model the hydrologic cycle of the Wabash River Basin. SWAT is a semi-distributed watershed modeling program developed by the Agricultural Research Service (ARS) of the U.S. Department of Agriculture [20,21]. SWAT model construction requires inputs of hydrography, topography, soils, and land cover. For this study, model construction was facilitated by the program ArcSWAT [22], a SWAT interface for the geographic information systems (GIS) software ArcGIS. Further model setup included the choice of the curve number method [23] for estimating surface runoff and the Penman-Monteith method [24,25] for estimating potential evapotranspiration (PET) and actual evapotranspiration (AET).
The following sections outline the main steps in setting up the Wabash Basin SWAT model, including the gathering and pre-processing the required model inputs (Section 2.2.1) and the definition of hydrological response units (Section 2.2.2). Additional details on model setup, including the incorporation of tile drains [26][27][28][29][30][31], the aggregation of pond and lake area by subbasin, and the parameterization of reservoirs, are included in supplemental Text S1, Figures S1 and S2, and Tables S1 and S2.

Model Inputs
SWAT input datasets, including topography, delineated subbasins, soil properties, and land use/landcover (LULC) were compiled from multiple governmental agencies. Basin topography, in the form of a 10 m (1/3 arc-second) digital elevation model (DEM), was obtained from the National Elevation Dataset (NED) via the National Map Seamless Server (https://viewer.nationalmap.gov/basic/). Delineated subbasin polygons were obtained from the United States Geological Survey's National Hydrography Dataset (NHD), with pre-defined streams and subbasins in the SWAT model corresponding to NHD 12-digit Hydrologic Unit Code (HUC) basins ( Figure S1). Soil properties were obtained from the State Soil Geographic Database (STATSGO), which lumps soils into four hydrologic groups. Soils within the Wabash Basin are dominated by hydrologic groups B and C (Table S5).
Two land use/landcover (LULC) datasets were used to parameterize SWAT: (1) the 2001 National Land Cover Data (NLCD) [32] and (2) the National Agricultural Statistics Service (NASS) cropland data layers [33]. The NASS dataset contains detailed spatial data on agricultural crop types-information that is not available in the NLCD dataset. NASS, however, has missing data due to cloud cover and has less detailed classification for non-agricultural lands. The distribution of different crop types is potentially important for future SWAT model applications. For example, corn and soybeans have different optimal growth temperatures and will, therefore, respond differently to continued climate warming. Therefore, a combined NLCD and NASS LULC dataset was created for the Wabash Basin (Table S4; Figure S3). The primary LULC classes are soybean and corn, encompassing 26.1% and 25.3% of the area within the Wabash Basin. Details on the pre-processing steps for the LULC inputs can be found in supplemental Text S1.
Historical meteorological data and climate projections were obtained from the University of Notre Dame [34,35]. These gridded time series have a daily time step and include maximum temperature, minimum temperature, precipitation, and wind speed. The climate change projections have been statistically downscaled using the hybrid delta method (HD) for a subset of ten different global climate models (GCMs) from Phase 5 of the Coupled Model Intercomparison Project (CMIP5) under representative concentration pathways (RCPs) 4.5 and 8.5 (Table 1). RCP 4.5 represents a medium stabilization scenario; and RCP 8.5 represents a very high baseline emissions scenario [36]. The gridded time series have a 1/16° spatial resolution and were aggregated at the subbasin level for use in the SWAT model.
For the historical baseline, a single model simulation was completed using the gridded observation-based meteorological dataset . Three future periods, representing 30-year windows centered on the 2020s (2011-2040), 2050s (2041-2070), and 2080s (2071-2100), were used for the climate change scenario modeling. The 10-member GCM ensemble combined with the two emissions scenarios (RCPs 4.5 and 8.5) lead to a suite of 60 SWAT model runs for the climate change analysis (10 GCMs × 2 RCPs × three future periods), plus one SWAT model run for the historical baseline. Model runs were completed using Indiana University's high-performance computing resources, with post-processing and model result visualization and analysis supported by Extreme Science and Engineering Discovery Environment (XSEDE) [37]. While the climate change projections are keyed to 30-year time periods, they have the same time series length as the baseline historical dataset due to the HD downscaling method, as described in [34,38]. For all model runs, the first 14 calendar years were used as the model spin-up period, resulting in 85 years (1929-2013) of daily model output. These relatively long time series have key advantages for estimating changes in extremes [38]; therefore, this dataset is well-suited for investigating changes in drought frequency and severity.
The 10-member GCM ensemble projects increases in temperature across all seasons, reaching a 5 °C increase in the mean annual temperature by the 2080s relative to the 1980s baseline under RCP 8.5. A comparison of annual and seasonal precipitation between the historical baseline and the three future periods, 2020s, 2050s, and 2080s, is presented in the results section, alongside the SWAT model water balance outputs. Further details on the projected changes in climate, the HD downscaling method, and the rationale for the selection of the 10-member GCM subset can be found in [34].

Hydrological Response Unit (HRU) Definition
In SWAT, hydrological response units (HRUs) are defined as unique combinations of land cover, soil, and/or slope classes within a subbasin. The HRU method is an effective way to simplify representation and simulation of watershed processes [39]. ArcSWAT allows users to specify two types of thresholds to define HRUs: percent-based and area-based. Higher thresholds result in fewer HRUs and thus shorter model run times. However, higher thresholds also result in greater generalization and thus a greater loss of information. Slope within the Wabash Basin is strongly correlated with landuse/landcover and soils; therefore, only one slope class was used, with slope set to the subbasin average. Based on an analysis of the trade-off between model simplicity versus information loss (see supplemental Text S1, Figure S4), area-based thresholds values were set as follows: land cover = 250 ha and soils = 650 ha, leading to a total of 6852 HRUs.

Model Calibration and Validation
After the initial data assimilation and model construction, the Wabash Basin SWAT model was calibrated and validated. Model calibration was completed in two stages. Stage 1 consisted of automated parameter regionalization, where a cascading methodology like the method outlined in [40] was used to set the calibration parameter values to an initial "best guess". Parameter regionalization was completed in nine subbasin sets, where upstream reaches were parameterized with the "best guess" parameter sets before moving to the next downstream gauges ( Figure S5, Table  S6). The goal of this regionalization was to account for regional variations in the performance of initial parameter values, especially with regard to the most sensitive parameters. For each round in this cascading scheme, 500 model runs were completed. Stage 2 of the model calibration consisted of a multi-site, multi-criteria sequential uncertainty fitting (SUFI) procedure, similar to the method outlined by [41]. For Stage 2, two rounds were completed with 850 model iterations each. This twostage calibration routine was completed using the R statistical programming software [42].
Both stages of the calibration routine use the same observed stream flow data and the same multi-criteria objection function, which are described in the following Sections 2.3.1 and 2.3.2. Additional details on the calibration methodology, including (1) parameter constraints for deep aquifer recharge [43] and ground water recession constants [44,45], (2) initial parameter ranges, (3) parameter modification with Latin hypercube sampling [46], and (4) baseflow separation [47,48] and the detailed use of the multi-criteria objective function, are included in Text S2, Figures S4-S6, and Tables S7 and S8.

Stream Flow Calibration Sites
Observed stream flow records for all gauging stations within the Wabash River Basin were downloaded from the United States Geological Survey (USGS) using the R package waterData [49]. A time series plot of total annual stream flow at the downstream main-stem gauge (ID 03377500) was used to select the calibration and validation periods for the model, with the goal of including high flow and low flow years in both periods and maximizing the number of active gauging stations. The 8-year period from 1993-2000 was chosen as the calibration period, with 2001-2012 serving as the validation period, and 1981-1992 serving as the warm-up period ( Figure S6).
As part of the initial data screening, a subset of gauging stations within the Wabash Basin was chosen based on: (1) record completeness-no missing data 1993-2012, (2) spatial distribution-avoid clusters of gauging stations, and (3) proximity to lakes reservoirs-gauging stations with reservoir and lake effects were removed. The final selection included 25 gauging stations located throughout the Wabash River Basin, with drainage areas ranging from 1383 km 2 to 74,164 km 2 ( Figure 2, Table  2).

Multi-Criteria Objective Function
The Nash-Sutcliffe efficiency criterion (NSE) [50] is often used to evaluate the performance of hydrological models. However, several shortcomings in the use of NSE as a single objective function have been pointed out [51] with one of the main problems being a greater emphasis on high flows. Therefore, other error metrics, including percent bias (PBIAS), root mean square error standard deviation ratio (RSR), and the volumetric efficiency criterion (VE) proposed by [51] were incorporated into a multi-criteria objective function (Table 3). VE was included because it represents the fraction of water delivered at the proper time. VE values range from 0 to 1, where 1 represents a perfect match. RSR and NSE are based on the minimization of the square of the residuals, while PBIAS and VE are based on the minimization of the absolute (VE) and relative (PBIAS) differences. These two error function types are complementary [52], and, thus, used together, help constrain parameter ranges, reduce model uncertainty, and increase model skill.
The "best" parameter set was chosen based on the minimization of a multi-criteria objective function, which incorporates the four different error metrics in Table 3. See Text S2 for additional details. For Stage 1 of model calibration, the objective function was calculated independently for each gauging station. For Stage 2, the objective function was calculated jointly for all 25 gauging stations by calculating a weighted mean for each error statistic, where the weights for each gauging station were equal to the ratio of the gauging station drainage (1383 km 2 to 74,164 km 2 ) area to the total drainage area upstream of all gauging stations (89,198 km 2 ).  1 The calibration target for error statistics RSR, PBIAS, and NSE are based on the threshold for "satisfactory" model performance at a monthly time step according to Table 4 in [53]. A daily time step is used here; thus, model simulations that meet the calibration target can be considered "good". VE is not listed in [53]; it has the same range as NSE and was therefore set to the same calibration target as NSE.

Quantification of Drought
To understand how climate change will impact drought, drought must first be quantified. A group of related drought indices-the standardized precipitation index (SPI) [54], the standardized soil moisture index (SSI) [55], and the standardized runoff index (SRI) [56]-are flexible, multi-scale indices that can be calculated over a range of timescales, typically 1 to 48 months. The multi-scale nature of these indices enables the quantification of both short-term and long-term drought characteristics. SPI, SSI, and SRI provide an assessment of meteorological, soil moisture, and hydrological drought, respectively. Hydrological drought is associated with deficiencies in the both surface water and groundwater; it is often diagnosed by streamflow drought. For this study, the quantification of hydrological drought focuses only on streamflow and does not consider groundwater storage.
The majority (66%) of the Wabash River Basin is cultivated agricultural land, and, in general, agricultural regions are the most sensitive to short-term droughts that coincide with critical crop development periods in July to August [19]. Therefore, SPI was calculated at 1-to 6-month timescales, and SSI and SRI were calculated at the 1-month timescale. The three drought indices-SPI, SSI, and SRI-were calculated using the same basic methodology: 1. The time series of watershed-averaged precipitation, soil water content, and runoff were obtained from the SWAT model output for the historical baseline model simulation.
2. The probability distributions of monthly precipitation, soil water content, and runoff were calculated separately for each month. Distribution fits were tested with the Shapiro-Wilk test.
The gamma distribution produced a satisfactory fit for precipitation. No satisfactory distribution was found for soil moisture or runoff. Therefore, percentiles for soil moisture and runoff were estimated empirically, following recommendations in [56]. 3. Monthly time series of cumulative probabilities were calculated using the fitted distributions (for baseline and climate change scenarios) and then converted to z-values, representing the number of standard deviations below or above the mean of a normal distribution with a mean of 0 and a variance of 1.
The standardized index values provide a measure of drought severity, with more negative values indicating more severe drought conditions. Drought frequency was calculated as the fraction of years with moderately dry (moderate drought) to extremely dry (extreme drought) conditions, i.e., standardized index value of less than −1 ( Table 4). The significance of changes in the frequency and severity of drought was then tested using the chi-square, Wilcoxon signed-rank, and Mann-Whitney U tests. The chi-square test was used to test for significant changes in drought frequency. The Wilcoxon signed-rank test [57] was used to test for significant shifts in the drought indices toward higher (wetter conditions) or lower (drier conditions) values. The Mann-Whitney U test [58] was used to test for significant changes in median drought severity.
The Wilcoxon signed-rank test is a nonparametric statistical hypothesis test to determine if there is a significant difference between two matched samples, where the matched samples correspond to the drought indices calculated from (1) the historical baseline and (2) the climate change scenarios. The Mann-Whitney U test is an alternate version of the Wilcoxon signed-rank test and is used to determine if there is a significant difference between to independent (i.e., unmatched) samples. For this study, the independent samples corresponded to the subset of drought index values (i.e., SPI, SSI, SRI values less than −1) for (1) drought events in the historical baseline and (2) drought events in the climate change scenarios.
To further investigate climate change impacts on the drought regime, the relationships between the different drought types were analyzed with bivariate correlation analysis and multiple linear regression. Identifying relationships between different drought types provides a basis for understanding how drought propagates through the hydrologic system, specifically the propagation of meteorological drought into soil moisture or hydrological drought. Therefore, SSI and SRI at the 1-month timescale were used as the dependent variables (predictands) with SPI at 1-to 6-month timescales used as the independent variables (predictors). Combinations of predictor variables were tested using automatic stepwise (combined forward and backward) regression analysis, with final model selection based on the Akaike Information Criterion (AIC) [59]. Extremely Wet

Results
The following sections present the model calibration and validation (Section 3.1) and the results of the climate change scenario modeling, including climate change impacts on the annual and intraannual water balance (Section 3.2) and the frequency and severity of drought (Section 3.3). Further analysis using multiple linear regression models to investigate how climate warming impacts the relationship between meteorological and soil moisture and hydrological drought is included in Section 3.4.

Model Calibration and Validation
The SWAT model was calibrated with the 1993-2000 period and validated with the 2001-2012 period. The error statistics for the calibration and validation periods for the "best" model simulation are included in Tables S9 and S10. For the calibration period, PBIAS was below the calibration target of ±25% for all 25 gauging stations, and the majority (23/25) of stations met the calibration goal for NSE and RSR. All calibration targets were met at the downstream main-stem gauging station (Wabash River at Mt Carmel-USGS gauge 03377500) in both the calibration and validation periods (Tables S9 and S10; Figures S8 and S9).
Overall, model performance in the validation period is similar to the calibration period. The main-stem gauging station (03377500) had a higher PBIAS value compared to the calibration period (−3.3% calibration, +7.0% validation). This pattern is similar for the other stream flow gauges; however, only three out of the 25 gauging stations have PBIAS values higher than the calibration target of ±25%. This tendency toward higher PBIAS values indicates an over-prediction of stream flow in the validation period compared to the calibration period. While the cause of the overprediction in the validation period is unclear, it may be due to changing land and water use patterns, which were not incorporated into the model. For this study, the model outputs for the full 61-member model ensemble are analyzed in terms of the sensitivity to climate change, focusing on the relative change rather than on absolute change. With this focus, the implications of the higher PBIAS values are minimal.

Climate Change Impacts on Annual and Intra-Annual Water Balance
A water balance analysis was completed for each of the 61 scenarios (one historical baseline; 60 climate change). The water balance components of interest for this study include the input precipitation and the model outputs of the amount of precipitation falling as snow, AET, soil moisture, ground water recharge, and stream flow. Water balance results for individual GCMs are not shown, but rather lumped by RCP (4.5, and 8.5) and period (2020s, 2050s, and 2080s) and reported by the inter-model spread and the ensemble mean for simplicity. Results are presented as the percent change from the historical baseline (((futurebaseline)/baseline) × 100) at both the annual ( Figure 3) and monthly (Figure 4) time scale.
Based on statistically downscaled climate projections, there is high certainty that regional temperatures will increase in the future [34]; however, the direction and magnitude of precipitation change is less certain. The inter-model variability of projected precipitation is high and increases from near future (2020s) to far future (2080s) (Figure 3a). In general, precipitation in winter (December, January, and February) is projected to increase while precipitation in late summer and fall (August, September, and October) is likely to exhibit no change or decrease slightly (Figure 4a).
Water balance variables that are closely linked to temperature, i.e., snow and AET, exhibit consistent changes and a low amount of inter-model variability (Figures 3 and 4). The amount of precipitation falling as snow decreases (Figure 3b) while AET increases (Figure 3c). As expected, the largest increases in AET occur in the far future (2080s) under the high emissions scenario (RCP 8.5). Seasonally, AET increases the most in winter (December to February) and decreases slightly or exhibits no substantial change during late summer and early fall (July to August; Figure 4b). The seasonal pattern in the AET change is related to seasonal water availability. Increasing temperatures cause PET to increase in all months (results not shown). AET, however, is dependent on both the PET and water availability. If little or no water is available, AET will be small or zero even when PET is high.
Soil moisture, ground water recharge, and stream flow are more closely linked to precipitation than to temperature. The mean annual values of soil moisture, ground water recharge, and stream flow exhibit both increases and decreases relative to the historical baseline, due to the high intermodel variability in precipitation (Figure 3d-f). Seasonally, ground water recharge, stream flow, and soil moisture exhibit a pattern of decreased water availability in summer (June, July, August) and fall (September, October, November; Figure 4c-e).  Table S11 for the tabular version.

Projected Changes in Meteorological, Soil Moisture, and Hydrological Drought
Soil moisture and hydrological droughts become significantly more frequent in the far future under both emissions scenarios, particularly during the summer months (Figures 5 and S10). Meteorological droughts, however, exhibit either decreasing frequency or no significant change, with relatively few GCMs showing a significant increase in meteorological drought frequency. Changes in monthly drought indexes, SPI, SSI, and SRI, mirror the patterns shown in the results for the monthly water balance presented in Section 3.2, with a shift toward drier conditions (lower SSI and SRI values) in summer and fall (Figure 6 and S11) and higher precipitation (higher SPI values) in winter.
Despite the increased frequency of soil moisture and hydrological drought ( Figure 5) and the general shift toward lower SSI and SRI values ( Figure 6) there is, in general, no significant change in the mean severity of moderate to severe drought events (Figures S12 and S13). Therefore, while soil moisture and hydrological droughts will likely be more frequent in the future, the average drought severity will not be significantly different from the historical baseline.  Figure S10 shows results for RCP 4.5.  Figure S11 shows results for RCP 4.5.

Climate Warming Impacts on Relationship between Drought Types
Soil moisture and runoff over 1-month timescales, represented by SSI-1 and SRI-1, exhibit the strongest correlations with precipitation over 3-to 4-month timescales (SPI-3, SPI-4; Table 5). In the historical baseline period, the magnitude of SPI, SSI and SRI are also closely linked, and near-normal precipitation results in near-normal soil moisture and runoff. However, as the climate warms, this relationship changes, and there is a shift toward lower soil moisture and lower runoff for the same near-normal amount of precipitation (Figures S14 and S15). This shift toward lower runoff and lower soil moisture indicates an amplification of soil moisture and hydrological drought due to climate warming. To quantify the degree of drought amplification, multiple linear regression (MLR) models were used. Since the relationships between SPI and SSI/SRI vary seasonally (Figures S14 and S15), separate MLR models were created for each month, for a total of 732 MLR models each for SSI and SRI (61 SWAT simulations × 12 months). The final MLR models account for an average of 71% and 78% of the variability in SSI and SRI, respectively ( Figure S16).
From these MLR models, the y-intercept represents the predicted SSI or SRI from an SPI of 0 and thus provides a simple estimate of the degree of drought amplification for each of the climate change scenarios. As expected, the y-intercept is non-significant and near zero for the baseline historical 1980s, when near-normal precipitation (SPI of 0) corresponds to near-normal soil moisture and runoff (SSI and SRI of 0). For the three future periods (2020s, 2050s, and 2080s), however, the y-intercept has a negative value and is significant at the 5% level in most months ( Figure S17), indicating that the shift toward lower soil moisture and lower runoff for the same near-normal amount of precipitation is statistically significant.
Plotting the y-intercepts from the MLR models versus change in monthly mean temperature shows that, in these model simulations, climate warming is driving the amplification of soil moisture and hydrological drought. Further, the rate of this amplification, in terms of the change in standardized drought index per °C (SI °C −1 ), can be quantified using the simple linear model fits shown on Figures 7 and 8. For soil moisture, drought amplification from climate warming is greatest in the months of February, March, and April. For runoff, it is greatest in the months of March and April. The implications of this are clarified by framing the amplification in terms of drought severity. In the historical baseline period, moderate meteorological droughts propagate into moderate soil moisture and hydrological droughts. A slope of −0.1 SI °C −1 indicates that, with 5 °C of warming, moderate meteorological droughts will be amplified into severe soil moisture and hydrological droughts. Further, under 5 °C of warming and a slope of −0.2 SI °C −1 , moderate soil moisture and hydrological droughts are likely to occur even with near-normal precipitation.

Discussion
While several recent studies have used hydrological models to investigate climate change impacts on meteorological, soil moisture, and hydrological drought [13,14,61], no previous studies have specifically focused on short-term drought events and identified intra-annual patterns in the drought response to climate change. In this study, a large ensemble of 61 hydrological models was combined with an analysis of monthly drought indices and MLR models. Climate change projections show increases in temperature and in the seasonality of precipitation, leading to increased frequency of soil moisture drought and hydrological drought. Hydrological modeling results show that climate warming amplifies soil moisture and hydrological drought, with the same amount of precipitation yielding significantly lower soil moisture and significantly lower runoff in the future than in the past.
The role of above-average temperatures in the amplification of soil moisture and hydrological drought has been documented by previous studies in California [16], central and western Europe [15], and the Iberian Peninsula [62]. In this study, we focus on a historically water-rich agricultural region within the US Midwest and quantify drought amplification in terms of the increase in drought severity per degree increase in monthly mean temperature. Results show that the degree of drought amplification varies seasonally. For soil moisture, the seasonal variation is related to the local hydroclimatology. With the same near-normal precipitation, climate warming amplifies soil moisture drought more in energy-limited months (e.g., March-April, Figure 7) when the PET and AET are near equal than in moisture-limited months (e.g., August-October; Figure 7) when the PET is substantially higher than the AET.
Climate warming also amplifies hydrological droughts most in spring (March-April, Figure 8). However, the impact of climate warming on hydrological drought is lowest in winter (December-January, Figure 8), likely due to the added complexity of snow accumulation and melt. Even with the seasonal variations, results clearly show that temperature plays a significant role in the amplification of soil moisture and hydrological droughts. Under 5 °C of warming, moderate meteorological droughts will be amplified into severe soil moisture and streamflow droughts for most months of the year. While the results presented here are region-specific and dependent on the local hydroclimatology, the methodology could be transferred to other regions and used to understand how the relationship between meteorological drought and soil moisture/hydrological drought will change under continued climate warming.
Like all hydrological models, the calibrated SWAT model used in this study is a simplified numerical representation of the natural flow system, and it cannot duplicate the natural flow system exactly. Uncertainties are inherent in model structure, input data, parameterization, and GCM selection. However, the model is physics-based and assessment with multiple error statistics indicates good performance (Tables S9 and S10). Attributing model uncertainty among the different sources is difficult and beyond the scope of this paper. Additionally, land use change and climate change mitigation were assumed to be constant in both the baseline and future climate change scenarios. Incorporation of land use change and climate change mitigation is beyond the scope of this paper but is a promising direction for future research.
Additional limitations of this study stem from aggregation of the water balance and drought indices to the watershed scale. The response of soil moisture drought and hydrological drought to climate change may vary spatially due to heterogeneity of the catchment's physical properties, e.g., land use, soil type, and geology. Additionally, this study only analyzed total soil moisture. Shallow soil moisture and deep soil moisture may exhibit different responses to climate change [63,64]. A useful extension of this work would be to (1) investigate if the physical properties of the catchment control the degree of soil moisture and hydrological drought amplification under climate warming and (2) determine if shallow and deep soil moisture exhibit different responses to climate warming.

Conclusions
In climate change projections, the uncertainty in precipitation is inherently higher than the uncertainty in temperature. Since precipitation dominates the hydrologic regime in water-rich regions, the uncertainty in projected precipitation carries through to uncertainty in the direction and magnitude of future change in the hydrological regime. The analysis presented here, however, separates out the role of temperature and shows that increased temperatures cause an amplification of soil moisture and hydrological drought. In general, for this study region, the frequency and severity of meteorological drought does not change significantly from the frequency and severity of meteorological droughts in the historical period. However, our analysis shows that while moderate meteorological droughts propagate into moderate soil moisture and hydrological droughts in the baseline period, this relationship changes in the future periods, and moderate meteorological droughts will propagate into severe soil moisture and hydrological droughts under 5 °C of warming.