Combined effects of predicted climate and land use changes on future hydrological droughts in the Luanhe River basin, China

Both climate and land-use changes can influence drought in different ways. Thus, to predict future drought conditions, hydrological simulations, as an ideal means, can be used to account for both projected climate change and projected land-use change. In this study, projected climate and land-use changes were integrated with the Soil and Water Assessment Tool (SWAT) model to estimate the combined impact of climate and land-use projections on hydrological droughts in the Lutheran River basin. We showed that the measured runoff and the remote sensing inversion of soil water content were simultaneously used to validate the model to ensure the reliability of the model parameters. Following calibration and validation, the SWAT model was forced with downscaled precipitation and temperature outputs from a suite of nine global climate models (GCMs) based on CMIP5, corresponding to three different representative concentration pathways (RCP2.6, RCP4.5 and RCP8.5) for three distinct time periods: 2011–2040, 2041–2070 and 2071–2100, referred to as early century, mid-century and late-century, respectively, and the land use predicted by the CA–Markov model in the same future periods. Hydrological droughts were quantified using the standardized runoff index (SRI). Compared to the baseline scenario (1961–1990), mild drought occurred more frequently during the next three periods (except for the 2080s under the RCP2.6 emissions scenario). Under the RCP8.5 emissions scenario, the probability of severe drought or above occurring in the 2080s increased, the duration was prolonged, and the severity increased. Under the RCP2.6 scenario, the upper central region of the Luanhe River in the 2020s and upper reaches of the Luanhe River in the 2080s were more likely to experience extreme drought events. Under the RCP8.5 scenario, the middle and lower Luanhe River in the 2080s was more likely to experience these conditions.


Introduction
Accelerated climate change can affect water availability globally. According to the International Panel on Climate Change (IPCC), extreme meteorological and hydrological events (e.g., droughts and floods) are expected to occur increasingly more frequently, which could cause more uncertainties and risks in river basins worldwide in the future (IPCC 2007(IPCC , 2014Fang et al. 2018a, b). In addition to climate change, alterations in land uses driven by anthropogenic activities lead to changes in water availability in variable ways (Nunes et al. 2011;Li et al. 2014;Wang et al. 2018). Hence, many studies have shown the combined impact of climate changes and human activities on the availability of water resources (Han et al. 2018;Kim et al. 2013;Ceola et al. 2014;Tong et al. 2012). Most of these studies adopted calibrated hydrological models forced with different climate and land-use projections to assess hydrological responses under the changing environment. These hydrological models generally vary from relatively simple lumped models (Liu et al. 2000) to modern process-based distributed models, such as the Soil Water and Assessment Tool (SWAT) model (Shrestha et al. 2017;Serpa 2015).
The approach most commonly used to evaluate climate change is through the application of global circulation models (GCMs), which are established to forecast future climate characteristics with different emission levels. However, GCM projections have low resolution due to the constraints of mathematical representations of atmospheric dynamics. To capture climate variability in regional hydrological simulations, projected GCM outputs processed by variable downscaling approaches were used to force process-based hydrologic models (Chattopadhyay and Edwards 2016;Stewart et al. 2015;Daggupati et al. 2016;Uniyal et al. 2015), in which statistical downscaling approaches were commonly used, such as the delta-change method (Tong et al. 2012;Dunn et al. 2012), weather generators (Seung-Hwan et al., 2013), bias correction (Felze, 2012) or regression-based downscaling. Different downscaling technique selections can lead to varying or even conflicting trends when the processed data are used to force hydrological models (Burger et al. 2012(Burger et al. , 2013. Therefore, it is important to choose an appropriate method to characterize future climate scenarios. It has been reported that decreasing vegetation could cause the annual mean discharge to increase in the context of intensive human activities (Andréassian 2004;Brown et al. 2005). Some researchers (e.g., Sala 2000) even believe that the influences of land-use alterations on hydrological responses could outweigh those caused by climate change. Current land-use projection methods generally contain generalized assumptions about future conversions (Dunn et al. 2012) and use modeling approaches based on socioeconomic and geographical driving factors (Kim et al. 2013). Land-use models are expected to predict the detailed land-use allocation driven by both geographic and socioeconomic factors. The variability between different land-use scenarios would profoundly influence the projections of hydrological responses (Piras et al. 2014;Stigter et al. 2014).
Previous studies on the influences of climate and land-use changes have mainly focused on water availability, but few studies have further investigated the evolutionary characteristics of drought. Drought is commonly considered a complex, multifaceted environmental hazard caused by climate anomalies. It usually arises from precipitation deficits, develops in hydrological systems and eventually threatens the socioeconomic system. Generally, drought is widely accepted to be divided into meteorological, agricultural, hydrological and socioeconomic droughts. Among these four types of drought, hydrological drought is most highly stressing due to its direct correlation with inadequate surface and subsurface water resources. In recent years, drought indices have developed into a popular means for drought detection due to the improvement of computational efficiency. To reflect the evolutionary characteristics of future drought under the changing environment, the selection of drought indices should consider both the ability to quantify drought conditions and the capability to depict drought propagation patterns under the same evaluation system. In this study, the standardized runoff index (SRI) (Zhang et al. 2015) was adopted to depict hydrological droughts.
The major objective of the present study was to evaluate the potential impacts of climate change integrated with land-use alteration on the evolution of future hydrological drought in the Luanhe River basin of China. The subobjectives included (i) making future projections of rainfall and minimum and maximum temperatures under a changing environment by the statistical downscaling model (SDSM); (ii) investigating future land-use changes in each catchment under a changing environment by using the CA-Markov model; and (iii) quantifying hydrological droughts under a changing environment by using the SRI. The findings can inform policy makers and resource managers in taking positive measures to cope with possible drought in the context of climate and land-use changes.

3 2 Study sites and data
The study area extends from 115°30′E to 119°15′E and 39°10′N to 42°30′ N. The total area of the Luanhe River basin is 44,600 km 2 (Fig. 1). This basin features a typical temperate continental climate that prevails with multiyear average temperatures varying from − 0.3 to 11 °C and average annual rainfall ranging from 400 to 700 mm/year. The intra-annual rainfall distribution is quite uneven, with more than 70% of the rainfall concentrated in summer (June, July and August). The occurrence of consecutive drought is more frequent due to decreasing runoff generation in the area and is becoming an enormous retardant to socioeconomic development. Major droughts in the Luanhe River basin have occurred in 1961, 1963, 1968, 1972, 1980 -1984, 2000, 2007 and 2009. In recent years, with the increasing demand for water resources due to global warming and social and economic development, the problem of water shortages and drought has become increasingly serious. Therefore, it is significant to assess future drought trends in a rapidly changing environment.
For downscaling purposes, reanalysis grid point data (2.5° × 2.5°) from 1961 ~ 2005 provided by the National Centers for Environmental Prediction/Nation Center for Atmospheric Research (NCEP/NCAR) were used for the study. The multimodel ensemble mean of 9 GCMs for CMIP5, as shown in Table 1, was used in this study to drive the SDSM to project future climate behavior. The choice of the models was dependent on the availability of GCM simulations for four representative concentration pathway (RCP) scenarios, namely  the RCP2.6, RCP4.5 and RCP8.5 scenarios, in this study area. The NCEP and CMIP5 data were obtained from the NOAA Physical Sciences Laboratory (https:// psl. noaa. gov/ data/ gridd ed/ data. ncep. reana lysis. html) and ESGF Portal at CEDA (https:// esgf-index1. ceda. ac. uk/ search/ cmip5-ceda/), respectively.
The hydrometeorological data used to drive the SWAT model included meteorological data, hydrological data and grid data. Information on meteorological, rainfall, hydrological gauges and grids is given in Table 2. The observed runoff from 1961 to 2007 was used to calibrate the SWAT model, and the rest of the recorded data (2008)(2009)(2010)(2011) were used to validate the model. The geographical location of the Luanhe River basin and relevant hydrometeorological stations used in this study is shown in Fig. 1.
The monthly soil water content data (0.25° × 0.25°) of GLDAS (Global Land Data Assimilation Systems) Noah_2.0 covering 1961-2010 were used to validate the established SWAT model.

Methodology
The methodology of this study consisted of the projection of future climate, hydrological modeling and land-use change modeling (Fig. 2). The SDSM and bias correction methods were applied to project the future climate, the SWAT model was used to simulate the runoff series, and the CA-Markov model was used to project the future land-use change in the Luanhe River basin.

SWAT model
The SWAT model is a dynamic model at the watershed scale that has a relatively perfect physical mechanism. It can set up a scenario model flexibly and conveniently to simulate the hydrological response under different scenarios. The hydrological process module, the soil erosion module and the pollution load module are the three main submodels of the SWAT model. The simulation of the hydrological process can be divided into slope hydrological processes and channel hydrological processes. The former determines the input of water, sediment, nutrients and chemicals of the main channel in each subbasin. The latter controls the transport of water, sediment, nutrients, pesticides and other substances in the river network. In this study, the hydrological process submodule of the model was used to simulate runoff and soil moisture changes.

SDSM method
The SDSM is a hybrid method of a multivariate linear regression method and a stochastic weather generator. The quantitative statistical relationship F can be described as follows (Chen 2000): where Y is the local predictand, and x i (i = 1 to n) is the ith large-scale atmospheric predictor.
The SDSM determines the precipitation occurrence and the precipitation amount through large-scale atmospheric variables as follows (Wetterhall et al. 2005): where ω t is the conditional possibility of precipitation occurrence on the tth day, ̂( j) t is the normalized predictor, the regression parameter α j can be obtained using the least squares method (LSM) method, and α t−1 is the t−1 day regression parameter. The last term of the equation α t−1 ω t−1 is optional. ω t is commonly compared to a uniformly distributed random number r t (0 ≤ r t ≤ 1) to determine whether precipitation occurs. The precipitation can be expressed by a z-score as follows: where Z t is the z-score on the tth day, and β j is the regression parameter. β t−1 and Z t −1 are the regression parameter and the z-score on the (t−1)th day, respectively. ε is a random error term.
Finally, precipitation y t can be expressed as follows: where ϕ is the normal cumulative distribution function, and F is the empirical distribution function of y t .

Quantile mapping (QM) of precipitation
The QM method is a nonparametric bias correction method that can be used without any assumption of precipitation distribution and has been widely used in previous studies  (Chen et al. 2013;Wilcke et al. 2013). It has been proven effective in correcting bias in the mean, standard deviation and wet-day frequency and quantiles. Equation (5) shows the precipitation adjustment expressed in terms of the empirical CDF (ecdf) and its inverse (ecdf −1 ) using QM as follows:

Distribution mapping (DM) of temperature
The idea of DM is to correct the distribution function of the raw data to agree with the observed distribution function. This can be done by creating a transfer function to shift the occurrence distributions of temperature (Sennikovs and Bethers 2009). It is used to adjust the mean, standard deviation and quantiles. Furthermore, it preserves the extremes (Themeßl et al. 2012). For temperature, the Gaussian distribution with a mean µ and a standard deviation σ is usually assumed to fit best (Schoenau and Kehrig 1990;Teutschbein and Seibert 2012): Similarly, the corrected temperature can be expressed as follows: where F N (⋅) and F −1 N (⋅) are the Gaussian CDF and its inverse, respectively; µ raw,m and µ obs,m are the fitted and observed means of the raw and observed precipitation series in a given month m, respectively; and σ raw,m and σ obs,m are the corresponding standard deviations, respectively.

CA-Markov method
CA-Markov is the combined cellular automata (CA)/Markov chain/multicriteria/multiobjective land allocation (MOLA) land cover prediction method. The Markov model mainly concentrates on quantifying land-use changes without detailed allocation of various land-use types in terms of the spatial extents (Wickramasuriya et al. 2009), whereas the CA model has a strong spatial conception. The CA-Markov model integrates the advantages of these two theories, which can yield a better effect in temporal and spatial land-use change simulations (Ji et al. 2009).
The specific procedures of using the CA-Markov model to simulate land-use changes are as follows: (1) calculating the land-use transition probability matrix to serve as the transformation rules adopted in the CA-Markov model; (2) determining the CA filters required to produce a clear sense of the space weighting factor; and (3) determining the starting point and the CA number of iterations.

Standardized runoff index
In this study, hydrological drought is described by the SRI, which is calculated similarly to the standardized precipitation index (SPI) but is instead applied to the runoff series. The calculation steps are as follows: (1) accumulate the runoff series at a given timescale; (2) fit a probability distribution to the specific runoff series; (3) estimate the cumulative density function (CDF) of an observed cumulative runoff volume; and (4) convert the cumulative probability to a standard normal function with zero mean and unit variance. For details about SRI computation, refer to the calculation of the SPI and SRI (Zarch et al. 2015;Ahmadalipour et al. 2017).

Performance evaluation criteria
In this study, the coefficient of determination (R 2 ), root mean square error (RMSE) and Nash-Sutcliffe coefficient (NSE) (Nash and Sutcliffe, 1970) were used to measure the performance of the SWAT and that of the SDSM (Fang et al. 2018a, b). R 2 is the square of the correlation coefficient and describes how much of the variance between the observed and simulated variables can be explained through a linear fit. The RMSE explains the difference between the observed and simulated variables. In other words, it describes the spread of error or the performance of the model. The RMSE is measured in the same unit used for the simulation or observed data. The NSE quantitatively describes the accuracy of the model output for hydrological variables. The NSE value can range between -1 and 1. These values are calculated as follows: where x obs,i is the ith observed precipitation or temperature (monthly), x sim,i is the ith simulated (raw GCM or downscaled) precipitation or temperature (monthly), and n is the number of data points. The quantitative accuracy test and the spatial accuracy test were used to measure the performance of the CA-Markov model. The quantitative accuracy and spatial accuracy can be defined as follows: where E 1 is the quantitative accuracy of class i land, and m iy and m ix are the simulated area and the current area of class i land, respectively. When E 1 is positive, it means that the simulated area of this type of land is larger than the current area; when E 1 is negative, it means that the simulated area of this type of land is smaller than the current area. The larger the absolute value of E 1 is, the higher the simulation accuracy is. E 2 is the spatial accuracy of class i land; c iy is the number of simulated correct raster cells; and c ix is the number of raster cells of class i land in the current land-use map. The larger the E 2 value is, the higher the spatial accuracy is.

Calibration and validation of the SWAT model
The model parameters of SWAT were calibrated and validated using the observed monthly runoff data from 8 hydrological stations (Fig. 1) from 1961 to 2011. Among them, 1961Among them, -1962 was the warm-up period used to reduce the influence of the initial conditions of the model at the initial stage of operation, 1963-1994 was the calibration period, and 1995-2011 was the validation period. The SUFI-2 method was used to analyze the uncertainty of the constructed SWAT model, and the R 2 , NSE and RMSE were used to evaluate the simulation accuracy and applicability of the model parameters. The calibrated parameter results of the SWAT model in eight subbasins of the Luanhe River basin are shown in Table 3. The observed and simulated values of eight hydrological stations were counted, and the corresponding relationship between runoff and rainfall is shown in Fig. 3. (Due to limited space, only the results of the CD, LY and LX stations are shown.) The calibration and validation results for the eight subbasins are shown in Table 4.
From the corresponding relationship between rainfall and runoff at eight hydrological stations (Fig. 3), the simulated runoff was higher in the months with more precipitation, indicating that the runoff was also greater when the precipitation was greater, and there was a positive correlation between them. In terms of the correspondence between precipitation and runoff, the spatial correspondence between weather stations and hydrological stations selected in the modeling process was relatively high. For both runoff quantity and runoff change trend, the simulated runoff series and the observed runoff series were well-fitted. The R 2 values for the calibration period and the validation period were all greater than 0.71 and 0.65, respectively, and the maximum values were 0.88 and 0.79, respectively, which were obtained at station LY. The NSE values were all greater than 0.68  1961 1966 1971 1976 1981 1986 1991 1996 2001 2006 2011 Precipitation ( 1961 1966 1971 1976 1981 1986 1991 1996 2001 2006 2011 Precipitation ( 1961 1966 1971 1976 1981 1986 1991 1996 2001 2006 2011 Precipitation (mm)

Precipitation
Observed streamflow Simulated streamflow Warm-up ValidaƟon and 0.63, respectively, and the maximum values were 0.85 and 0.79, respectively, which also appeared at LY station. The RMSEs were all less than 14.56 and 10.27, respectively, and the minimum values were 1.49 and 1.44, respectively, which also appeared at LY station (Table 4). The simulation effect of the peak flows in summer was poor, and the simulation value showed some underestimation, which may be related to the fact that the reservoir discharge in summer was not considered. The simulation effect for the validation period was slightly lower than that of the calibration period, which may be related to the different simulation accuracies of the SWAT model for wet and dry years, and the simulation accuracy in dry years was lower than that of wet years (Duan et al., 2014). Another possible reason was that the underlying surface of the basin changed during the verification period compared with the calibration period. Although the simulation effect in the verification period was lower than that in the calibration period, it met the requirements of model accuracy. The constructed SWAT model in this paper had a good simulation effect on the runoff process and presented good adaptability in the Luanhe River basin. It is feasible to simulate the hydrological process model in the Luanhe River basin by the SWAT model, which can be used for further research. By calculating the Spearman correlation coefficients for all subbasins between the soil moisture content of GLDAS Noah 2.0 and the soil moisture content simulated by the SWAT model, the simulation effect of the soil moisture content of the SWAT model was evaluated. The results of the correlation analysis are shown in Fig. 4. As shown in Fig. 4, the correlation coefficient between the assimilated soil moisture content and simulated soil moisture content in each subbasin was greater than 0.58, with a maximum correlation coefficient of 0.86, which appeared in subbasin #19. In general, the simulation accuracy upstream of the Luanhe River was significantly higher than that downstream. The results of the correlation analysis showed that the SWAT model could accurately depict the variation characteristics of the soil moisture content in the Luanhe River basin from 1961 to 2010.

Evaluation of the SDSM model for the calibration and validation periods
The 45-year observed data were used for the calibration  and validation (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) of the SDSM. The performance of the SDSM was evaluated by seven precipitation indices and two widely used temperature indices. These indices provided detailed 1 3 marginal information about the intensity, duration and frequency of precipitation series, rather than the only simple mean statistics.
The indices also provided detailed information about the intensity, duration and frequency of the time series. Table 5 gives the definition of the precipitation index used in this  Only the calibration and validation results of BLN station are presented graphically as a representative example because of the size of the paper. Figures 5, 6 present the monthly, seasonal and annual simulation effects of the seven precipitation indices for BLN station during the calibration and validation periods. It was obvious that the model showed good performance in the calibration and validation periods by means of the mean, variance, wet-day and Cday indices, and the Cwet; additionally, the 90th percentile and Rx5day indices were poorly matched between the observation and simulation data. The performance of the downscaling models was also assessed using standard statistical approaches, namely the R 2 , RMSE and NSE. Figure 7 shows the ensemble average of the R 2 , NSE and RMSE between the observed and simulated series evaluated daily for each station during the calibration and validation periods. It was obvious that the values of R 2 and NSE for the 45 subbasins were all larger than 0.35 for the calibration period, and they were all larger than 0.3 for the validation period. In addition, the RSME values were all smaller than 9 for the calibration and validation periods. The results showed that the simulation accuracy of the SDSM for precipitation was relatively low, and the reasons may be as follows. Compared with the influence mechanism of temperature and other meteorological phenomena, the stochastic processes affecting the formation of precipitation were more complex. Therefore, even though a possibly highest resolution GCM model was used as the input, the SDSM could not accurately simulate the amount and location of precipitation. The results of other research showed that the R 2 values of downscaling daily  (Wilby et al. 2002;Hessami et al. 2008), and the simulation results in the Luanhe River basin were close to these values. Therefore, the simulation ability of the built SDSM model to reproduce precipitation in this study was sufficient. However, the accuracy simulated in this study could not meet the requirements of the hydrological models. Therefore, bias correction was necessary before the output results of the SDSM model were used for hydrological simulations. In this study, the bias correction method of QM was performed on rainfall data downscaled using the SDSM model. The results of the bias correction of precipitation are presented in Fig. 8. As shown in Fig. 8, the correction methods adjusted the biases in the output precipitation of the SDSM model, and the corrected precipitation had quantile values similar to those of the observations. Figures 9, 10 present the mean and variance values of the maximum and minimum daily temperature at BLN station during the calibration and validation periods. As shown in Figs. 9, 10, the mean value and variation in the maximum and minimum daily temperatures can be accurately simulated by the SDSM for the calibration period. For the validation period, the mean of the maximum temperature was slightly underestimated, and the variation was slightly overestimated for all four seasons. The mean of the minimum temperature was slightly underestimated for summer and slightly overestimated for the other seasons, and the variation was slightly overestimated for all four seasons. Figures 11, 12 show the ensemble averages of the R 2 , NSE and RMSE between the observed and simulated maximum and minimum daily temperatures at each station during the calibration and validation periods. It was obvious that the values of R 2 and NSE for the 45 subbasins were all Fig. 6 Precipitation indices of observation and simulation for the validation periods at the representative station of BLN larger than 0.95, and the values of RSME were all smaller than 3, which indicated that the statistical downscaling procedure performed well for the maximum and minimum temperatures. To improve the simulated accuracy of the SDSM, the DM method was used for the maximum and minimum temperatures downscaled using the SDSM model. The corrected results are shown in Figs. 13,14. From Figs. 13,14, it was obvious that the temperature correction method adjusted the biases in the maximum and minimum temperatures and that the mean of the corrected temperature was closer to the values of the observations. The bias correction factors seemed to provide a good improvement in the correction of the maximum and minimum temperature from the SDSM model prior to their use in hydrological models.
The uncertainty of climate change impact assessments includes future emissions scenarios, climate models, downscaling technology, hydrological models, etc. Other studies by domestic and foreign scholars have found that compared with emissions scenarios and downscaling methods, the choice of GCM is the biggest factor leading to the uncertainty of climate change impact assessments. When evaluating the impact of climate change on future hydrology and water resources in a region, biased conclusions are likely to be drawn if only the scenarios simulated by one model are taken as the basis. Therefore, nine climate models commonly used in China were selected in this paper, and multiple climate models were collected using aggregation technology to reduce the uncertainty of a single climate model. Three scenarios representing low emissions (RCP2.6), medium emissions (RCP4.5) and high emissions (RCP8.5) were selected to reduce the uncertainty of carbon emissions scenario prediction. The downscaling results were corrected by the QM and DM methods, Fig. 7 Performance statistics of the downscaled daily precipitation from the SDSM model for the calibration and validation periods at 52 precipitation stations which were applicable to precipitation and temperature, respectively. The uncertainty of each link was reduced by multiple climate models, different carbon emissions scenarios and deviation correction methods to ensure the rationality and practicability of the climatic prediction results and hydrological simulation results in the Luanhe River basin.

Land-use change projection using the CA-Markov model
The transfer probability matrix of land-use change can not only quantitatively explain the conversion between land-use types but also reveal the transfer rate among different types. The transfer matrix of land-use area in the Luanhe River basin from 1990 to 2000 could be obtained by applying the CA-Markov model, as shown in Table 6. The starting point is the land use of 1990, and the number of iterations is 10. From 1990 to 2000, the major land-use types in the Luanhe River basin underwent complex mutual transformation. The areas of dry land, paddy fields, water bodies and unutilized land decreased (a) Calibration period (b) Validation period Fig. 8 Q-Q plots of the observed and simulated precipitation before and after bias correction for the calibration and validation periods at the representative station of BLN by 3.03%, 3.80%, 2.67% and 2.29%, respectively, and the areas of grassland, forest and construction land increased by 0.62%, 1.67% and 2.87%, respectively. Therefore, the dry land area decreased the most, and most of the land transferred out was forest and grassland, which was the result of returning farmland to forest and grassland under the guidance of government policies. Affected by the rapid development of economic construction and urbanization, the area of construction land increased rapidly, and the main sources of transfer were dry land and grassland. The water body areas decreased by 17 km 2 from 1990 to 2000. On the one hand, this change reflects that with the increase in urban expansion and agricultural intensification, the amount of surface water absorbed becomes larger, leading to lake and river reclamation. On the other hands, it is mainly affected by climate change. Wang et al. (2015) found that the Luanhe River basin has had a drying tendency over the past 50 years, especially since the 1980s. The drought index in the Luanhe River basin presented continuous negative values many times, and after 1999, this basin experienced more severe drought. This result was roughly consistent with the change trend of the water area in this paper, showing that the impact of climate change on the water area was significant (Yinglan et al., 2019).
The land use in 2005 and 2010 was projected using the CA-Markov model and compared to the actual land-use maps of 2005 and 2010 (Table 7). The simulated land-use maps of 2005 and 2010 were very similar to the observed land-use maps of 2005 and 2010, which reflected that the results simulated by the CA-Markov model had good reliability and applicability. The relative error of area and spatial distribution of each land use category for the calibration period were less than 10% and 20%, and for the Fig. 9 Mean and variance of the observed and simulated maximum temperatures for the calibration and validation periods at the representative station of DL validation period, the relative errors of the area and spatial distribution were less than 10% and 25%, respectively, showing the total accuracy of the land-use simulation to be reliable; hence, the CA-Markov model was used to simulate land-use changes in the Luanhe River basin in the 2020s, 2050s and 2080s.
On the premise of proving the feasibility and accuracy of the model, the land use in 2010 was taken as the initial state, and the transfer probability matrix of land use from 1900 to 2000 was adjusted based on the land-use types in 2005 and 2010. The land-use structure in the next 25 years, 55 years and 85 years was predicted by the model, and the results are shown in Fig. 15.
According to the predicted results, the changes in major land-use types in the Luanhe River basin in the future period (2020s, 2050s and 2080s) were basically consistent with the change trends in 1990-2000. Dry land, paddy fields, water bodies and unutilized land continued to decrease; the reduction rates of dry land, water bodies and unutilized land were faster; and the amounts of forest and grassland continued to grow slowly (Fig. 15). Compared with 2010, the forest area increased the most, followed by grassland in 2025, and forest accounted for 40.27%, 40.48 and 40.65% of the total area in 2025, 2055 and 2085, respectively. Dry land decreased the most, followed by unutilized land, and dry land accounted for 21.02%, 21.00 and 20.96% of the total area in 2025, 2055 and 2085, respectively. The largest annual increase was in construction land. In contrast, the largest annual decrease was in unutilized land. In addition, the change ranges of forest, grassland and paddy field area were small and basically stable.

Impacts of climate change and land-use changes on hydrological drought
The daily climate data obtained through the downscaling method and the land-use data obtained through the CA-Markov method in the future periods were substituted into the calibrated SWAT model to simulate the runoff series of each subbasin under future climate scenarios. On this basis, the SRI was calculated, and then, the characteristics of drought variation in the future periods were identified using run theory. Figure 16 shows the variation in the drought duration, drought severity and drought frequency in the Luanhe River basin during the baseline period and in future periods. In both periods, the higher the drought level was, the smaller the drought duration and drought frequency were, indicating that the duration of severe drought was short and the probability of occurrence was low. The duration of mild drought in the 2020s under the three emissions scenarios and in the 2050s under the RCP2.6 and RCP4.5 emissions scenarios was higher than that in the baseline period. The severity of mild drought in the 2020s under the RCP2.6 emissions scenario and in the 2020s and 2050s under the RCP4.5 emissions scenario was all higher than that in the historical baseline period. Under the (a) Calibration period (b) Validation period Fig. 11 Q-Q plots of the observed and simulated maximum temperatures before and after bias correction for the calibration and validation periods at the representative station of DL three emissions scenarios, mild drought occurred more frequently in the next three periods (except in the 2080s under RCP2.6) than in the historical baseline period, indicating that the probability of mild drought increased in the future. Under the RCP2.6 emissions scenario, the moderate drought duration in the 2020s and that under the RCP8.5 emissions scenario in the 2080s was higher than that in the baseline period. Under the RCP8.5 emissions scenario, the moderate drought severity in the 2080s was greater than that in the baseline period. Under the three emissions scenarios, the frequency of drought during the next three periods was less than that in the historical baseline period; that is, the probability of drought in the future period decreased. Under the RCP8.5 emissions scenario, the extreme drought duration of the 2080s was longer than that of the baseline period, and the drought severity and frequency increased. In general, the drought characteristics of light drought changed more greatly than that in the baseline period. Under the RCP8.5 emissions scenario, the probability of severe drought or above occurring in the 2080s increased, the duration was prolonged, and the severity increased. Figures 17, 18 show the spatial distribution of the drought duration, drought severity and drought frequency of extreme drought in the future under three emission scenarios for the Luanhe River basin. From Figs. 17 and 19, it can be seen that under the RCP2.6 and RCP4.5 scenarios, the duration of extreme drought in the 2050s was shorter, and the drought severity was smaller. Under the RCP2.6 and RCP8.5 scenarios, the duration of extreme drought was longer, and the drought severity was greater in the 2080s. The long drought duration and large drought severity during the historical baseline period Under the RCP2.6 scenario, the long drought duration and great drought severity in the 2020s, 2050s and 2080s were mainly distributed in the upper central region, the western region of the lower reaches, and the upstream region and the western region of the lower reaches of the Luanhe River, respectively. Under the RCP4.5 scenario, the long drought duration and great drought severity in the 2020s, 2050s and 2080s were mainly distributed in the middle and upper reaches, the western region of the middle reaches, and the western region of the lower reaches, respectively. Under the RCP8.5 scenario, the long drought duration and great drought severity in the 2020s, 2050s and 2080s were mainly distributed in the middle reaches, the upper central region and the western region of the lower reaches, and the middle reaches of the Luanhe River, respectively. Under the RCP2.6 scenario for the 2080s and the RCP4.5 scenario for the 2020s, the Luanhe River upstream was prone to extreme droughts, with a longer drought duration than that in the historical baseline period. Under RCP8.5 for the 2080s, the middle and lower reaches of the Luanhe River were prone to extreme droughts, which had longer drought durations and greater drought severities. It can be seen from Fig. 18 that under the RCP2.6 scenario in the 2020s, the upper central region of the Luanhe River was more likely to experience extreme drought events than during the historical baseline period, and in the 2080s, the upper reaches were more likely to experience these effects. Under the RCP8.5 scenario in the 2080s, the middle and lower reaches were more likely to experience extreme drought events.

Conclusions
This paper investigated the combined impacts of future climate change and land-use change on hydrological drought in the Luanhe River basin, China. To obtain reliable precipitation and the T max and T min sequences under various climate change scenarios, the SDSM model was used to downscale future climate variables, and the QM and DM methods were used to correct the output precipitation and temperature of the SDSM model. The CA-Markov model was applied to predict the land-use conditions in the 2020s, 2050s and 2080s. The future climate scenarios and land-use conditions were input into a well-validated SWAT model to simulate future hydrological processes. Based on the simulated runoff series, the SRI was used to characterize the hydrological drought in the 2020s, 2050s and 2080s under the RCP2.6, RCP4.5 and RCP8.5 scenarios. The most notable conclusions of this study were as follows: (1) The constructed SWAT model had a good simulation effect on the runoff process. The observed and simulated values were in good agreement, indicating that the model calibrated by a set of optimized parameters could be used to study the response of hydrological drought to climate and land-use change in the Luanhe River basin.
(2) The SDSM model will suffice to reveal the statistical relationships between the largescale atmospheric variables and the observed regional-scale meteorological factors, except for its less satisfactory performance for precipitation. The bias correction factors seemed to provide good improvement in the correction of the precipitation and maximum and minimum temperatures from the SDSM model. Hence, the results were acceptable for practical use.
(3) Land-use change analysis with the CA-Markov model projected that forest would increase the most, followed by grassland, and that dry land would decrease the most, followed by unutilized land. The largest annual increase was in construction land. In contrast, the largest annual decrease was in unutilized land, and the change ranges of forest, grassland and paddy fields were small and basically stable. (4) Under the three emissions scenarios, mild drought occurred more frequently during the next three periods (except in the 2080s under RCP2.6) than during the historical  baseline period, indicating that the probability of mild drought increased in the future period. Under the RCP8.5 emissions scenario, the probability of severe drought or above in the 2080s increased, the duration was prolonged, and the severity increased. Under the RCP2.6 scenario in the 2020s, the upper central region of the Luanhe River was more likely to experience extreme drought events, and in the 2080s, the upper reaches were more likely to experience these conditions. Under the RCP8.5 scenario in the 2080s, the middle and lower reaches were more likely to experience extreme drought events.