Evaluation of hydrological responses to climate change for a data-scarce mountainous watershed in Taiwan

Improved living standards have resulted in a substantial increase in the human population; overall, this has resulted in a greater demand for various resources (i.e., food, energy, and water). For the proper management of water resources, predictive hydrological models are necessary to quantitatively map their dynamics (Tsegaw et al. 2020). However, global climate change poses challenges to managing water resources. Liu et al. (2016) showed that global trends in regional annual precipitation were primarily driven by the changes due to the top 30% of heavy precipitation events during the major global warming period from 1979 to 2013. Sharma et al. (2018) noted that the discharge in small catchments is highly sensitive to temperature. Due to the relatively smaller spatial range of catchments, storms tend to encompass an entire catchment, resulting in complete saturation of the soil, thereby impacting streamflow. Therefore, changes in precipitation and temperature can significantly impact hydrological responses, which may lead to potential disasters due to water resource scarcity (IPCC 2022).

Predicting the responses of a hydrological system to climate change is important within the context of effectively managing water resources. However, challenges in understanding the hydrological responses of watersheds to climate change are often caused by the complex interactions of the system components, different scales of hydrological processes, difficulty in accurately measuring hydrological components, and insufficient long-term hydrological records. Owing to the complex interactions of the natural and anthropogenic influences in the hydrological cycle, it is difficult to isolate individual factors that shape it. Thus, mountainous watersheds and rural areas are ideal for the exploration of natural hydrological responses to climate change; this is primarily due to the fact that most of them have experienced minimal anthropogenic interferences. However, data on mountainous watersheds are usually scarce or not comprehensive (Afshar et al. 2017). Nonetheless, hydrological analyses have been conducted on data-scarce watersheds (Huang & Yang 2015; Liang et al. 2015; Tegegne et al. 2017).

Hydrological analyses can be classified into two groups, namely, physics-based (knowledge-driven) models and data-driven (empirical) models (Zheng et al. 2018). Data-driven models have been widely used owing to their simplicity and flexibility (Sivakumar & Berndtsson 2010). Specifically, data-driven models, such as the unit hydrograph method, various statistical models (i.e., linear regression, multi-linear, and ARIMA), and machine learning models (Anaraki et al. 2020; Farzin & Anaraki 2021; Farzin et al. 2022; Kadkhodazadeh et al. 2022), rely on a large amount of training data to describe the responses of the modeled system (Hurst 1951; Solomatine & Ostfeld 2008; Conoscenti et al. 2014; Taormina & Chau 2015; Pourghasemi et al. 2017; Bai et al. 2019). On the other hand, physics-based models are based on an understanding of system behavior (Zahabiyoun et al. 2013), and simulations of the hydrological responses are based on the first-order principles from the physics of a phenomenon or system (Liu et al. 2020). This type of model can be further subdivided, depending on the treatment of the system, into distributed or lumped models. Distributed models describe spatiotemporal hydrological responses, whereas lumped models provide representative hydrological characteristics for the entire watershed. Distributed models are often used for rigorous hydrology modeling. To date, several well-developed distributed physics-based models have been developed, including GR4J (Perrin et al. 2003), HSPF (Borah & Bera 2003), and the soil and water assessment tool (SWAT) (Bieger et al. 2017). The advantages of using distributed physics-based models include their excellent forecasting ability and the successful modeling of detailed spatial information using limited hydrological data (Balica et al. 2013). Therefore, physics-based models may be preferred for watersheds with scarce data.

Spectral analysis is a powerful technique for measuring the persistence of hydrological time series (Zhou et al. 2020). Several studies have characterized the fractal behavior of hydrological responses through spectral analysis (Bai et al. 2019). The fractal phenomenon, found across various fields (López-Ortega et al. 2020), was first reported by Hurst (1951) during an investigation of the discharge time series of the Nile River within the design framework of the Aswan High Dam. Kirchner et al. (2000) showed that the fractal scaling of the stream chemistry may range over three orders of magnitude, even when rainfall chemistry appears as a white noise spectrum. Chen & Hsu (2007) evaluated the influence of the fractal geometry of porous media on groundwater flow and solute transport. Habib (2020) implemented fractal analysis to determine the long-term hydrological memory features. Sun et al. (2019) conducted a spectral analysis to determine the fractal dimensions of a groundwater system. Fractal analysis of a time series is based on quantifying the degree of fluctuation around the overall trend of the data over time. Such an approach can measure the complexity of a system and be adopted in this study.

Climate change may result in significant changes in the responses of hydrological systems in various areas of the world, including Taiwan. Guo et al. (2018) presented the frequency of extreme hydrological events and surface and groundwater changes in Taiwan owing to climate change. Hsu et al. (2007) quantitatively evaluated the impact of climate change on groundwater resources in the southern plains of Taiwan. Yeh et al. (2015) used the Mann–Kendall test to analyze the dynamics of streamflow data from 12 gauges, revealing a crucial decrease in the streamflow in Taiwan due to climate change, which endangers its surface water supply. Wang et al. (2020) indicated a possible remote connection in the hydrological response of island-scale subtropical and continental-scale frigid zones.

Similar to other areas of the world, Taiwan also encounters crises related to uncertain water resources due to climate change. For example, from late 2020 to June 2021, Taiwan experienced a severe water shortage because no single typhoon hit Taiwan in the expected typhoon season in 2020. Consequently, the country's reservoirs could not be replenished for the first time in 56 years. To conserve water, Taiwan stopped irrigating 24% of the total planted rice crop areas in 2021 (Lin & Syngle 2021). In addition, Taiwan recorded its hottest month in May of 2021; specifically, this was the hottest month since 1947. Thus, exploring the impact of climate change on hydrology is important for managing routine and emergency water resources.

Socioeconomically, mountainous watersheds have been shown to be important within the context of providing water and hydroelectricity. However, the current understanding of the hydrologic response to climate change is limited, and this is largely due to a scarcity of the necessary data. Therefore, this study aims to (1) explore the impacts of climate change on watershed hydrology and (2) perform hydrological analyses using an integrated approach for data-scarce watersheds. To the best of our knowledge, this is the first study to characterize and predict the hydrological response to climate change in a mountainous watershed in Taiwan. These results will provide an understanding of the distribution of water components under climate change; such insights are valuable within the context of water resource management. In addition, the proposed approach can be applied to data-scarce watersheds around the world to measure the impact of climate change on future water scarcity. The remainder of the paper is structured as follows: the methodology section describes methods tied to the use of physics-based hydrological, groundwater, and general circulation models (GCMs). In addition, spectrum analysis is introduced, and the power spectral exponent is described; the study site section provides information on the Wuling Mountain Watershed. The hydrological data are described in the following section. In the Results and Discussion section, the integrated approach is applied to the Wuling Mountain Watershed. The results of the SWAT with modular three-dimensional groundwater flow (SWAT-MODFLOW) modeling, spectrum analysis, and hydrological simulations are presented. Power spectral exponents are calculated for hydrological components of the future scenario and compared to those of the baseline period. The impact of climate change on the hydrological response of the Wuling Mountain Watershed and its implications for future water resources are discussed. Additionally, the limitations and assumptions of the study are discussed. Finally, the study is concluded.

SWAT-MODFLOW simulates the surface water–groundwater interactions by integrating the SWAT with the modular three-dimensional groundwater flow model (MODFLOW) (Kim et al. 2008; Jafari et al. 2021). SWAT simulates the hydrological processes and solute transport above the aquifer, whereas the MODFLOW models the groundwater flow in an aquifer (Park et al. 2019). The recharge calculated by SWAT serves as the upper boundary of the MODFLOW for the simulation of the groundwater flow in aquifers.

The MODFLOW groundwater flow numerical model (McDonald & Harbaugh 1988) is one of the most widely used groundwater models. It can be used to simulate a multi-layer aquifer system in steady, transient, confined, or unconfined aquifers. Within the groundwater system, MODFLOW can simulate external disturbances, including pumping, regional replenishment, evaporative dispersion, water level changes, and regional drainage.

SWAT and MODFLOW are integrated as SWAT-MODFLOW. In SWAT, the smallest calculation HRU is a unique combination of land use/cover, soil, and slope and does not have a designated geographic location. Therefore, HRUs are disaggregated into individual polygons with specific geographic locations. These disaggregated HRUs are then intersected with the grid cells in MODFLOW to share common variables for SWAT and MODFLOW modeling. In MODFLOW, the river used to estimate the volumetric flow exchange rates between the aquifer and stream is intersected with the sub-basins in SWAT to transfer the groundwater return flow rates to the corresponding sub-basin stream.

A smaller RMSE indicates a closer difference between an observation and an estimation.

Moriasi et al. (2007) and Parajuli et al. (2009) provided the following criteria for model performance: very good performance: 0.75 < NSE ≤ 1 and 0.9 ≤ R² ≤ 1; good performance: 0.65 < NSE ≤ 0.75 and 0.75 ≤ R² ≤ 0.89; and satisfactory performance: 0.5 < NSE ≤ 0.675 and 0.5 ≤ R² ≤ 0.74. A lower RMSE value indicates better performance of the model. The SWAT-MODFLOW has been successfully applied to investigate the surface water–groundwater interaction of an irrigated alluvial fan in central Taiwan (Ke 2014).

The latest version of the statistically downscaled projections obtained from the Taiwan Climate Change Projection Information and Adaptation Knowledge Platform (TCCIP; https://tccip.ncdr.nat.gov.tw/) provides 5 km by 5 km grid-based projections for precipitation and temperature based on the CMIP5 (the fifth phase of the Coupled Model Intercomparison Project). The CMIP6 version of statistically downscaled projections is expected to be released in early 2023. The CMIP5 experiments were conducted using several GCMs under four representative concentration pathways (RCPs) of greenhouse gas emissions (i.e., RCP2.6, RCP4.5, RCP6.0, and RCP8.5) (IPCC 2013). The baseline period used to investigate the impact of climate change is from 1986 to 2005, and the future period is from 2021 to 2040 (WRA 2016). Based on the downscaled information, 200 years of precipitation and temperature data were generated under various climate change scenarios. The methods used for weather generation are as follows:

Precipitation generation is composed of the following steps: (1) simulate the occurrence of wet and dry days using the Markov chain process and (2) generate the rainfall amount for wet days using a specific probability distribution. In Step 1, a first-order, two-state Markov chain was applied to describe the transitions of a system state according to the transition probability matrix. The model assumes that there are only two possible system states (i.e., wet and dry), where the current system state depends on the last system state. Each element of the transition probability matrix can be derived based on the observed daily rainfall data. Due to seasonal features, the matrix for each month was calculated. For example, P_m(w) is the probability of a wet day in a month m and P_m(w|w) and P_m(w|d) are the conditional probabilities of a wet day when the previous days are wet and dry days, respectively, in month m. A uniform random variable x between 0 and 1 is employed to assign the state of the first day. When the variable x is less than or equal to P_m(w), the system state is defined as a wet day; otherwise, the system state is defined as a dry day. The system states of days are assigned based on a transition probability matrix. For example, if the state of the previous day was wet, the current day will be generated as a wet day when the random variable x is less than or equal to P_m(w|w); otherwise, the current day is simulated to be a dry day. In Step 2, after deciding the system state using the Markov chain, the rainfall amount was sequentially generated based on a specific probability distribution. Several probability distributions are applicable when generating the daily rainfall amount, including exponential (Tung & Haith 1995), Weibull (Yu & Polycarpou 2002), two-parameter gamma (Corte-Real et al. 1999), and mixed exponential (Woolhiser & Roldan 1986) distributions. Among these probability distributions, the Weibull distribution appropriately approximates the daily rainfall in Taiwan (Yu et al. 2002). It produces reasonable mean and standard deviation values for rainfall. Thus, a Weibull distribution was adopted in this study to generate the rainfall amount. The mean and variance of the observed daily rainfall data were calculated for each month to obtain the Weibull distribution parameter. Subsequently, a cumulative distribution function for the Weibull distribution was linked with a uniform random variable x to generate the amount of rainfall on a wet day.

The rainfall generation model can further process downscaled projections to generate rainfall data for various climate change scenarios. For example, the mean of the observed daily rainfall can be adjusted (e.g., −10%) according to the downscaled projections. However, the remaining rainfall characteristics (dry/wet status and variance) are assumed to be identical to the observed rainfall.

Equation (10) shows that the energy of the time-domain signal is equal to the energy of the frequency-domain transform, ⁠.

β describes the shape complex of an object and provides a statistical index of complexity, which compares the changes in patterns with their measurement scale (Malamud & Turcotte 1999). When β > 1, the time series is nonstationary with strong persistence. When 1 > β > 0, the time series is stationary and weakly persistent. When β = 0, the time series is uncorrelated. When β < 0, the time series is stationary and exhibits antipersistence behavior.

The Wuling hydrological station is located at a latitude of 24°08′14″N and a longitude of 121°16′46″E, at an altitude of 3,087 m. Hourly meteorological and hydrological data, including precipitation, temperature, relative humidity, and streamflow, were collected from 14 March 2013 to 22 July 2015. Two cluster wells (DHW-06) (121°16′54″, 24°08′35″) are at an elevation of 3,023 m. The daily groundwater levels were measured and recorded. The two wells are 3 m apart. The shallow well is 20 m below the ground surface (BGS) and the deep well is 100 m BGS. The deep and shallow wells have diameters of 2 and 6 inches, respectively. The screening depth of the well is 5–17 m BGS for the shallow well, while there are two screens located at 61–67 and 88–94 m BGS for the deep well.

The hydrological model was calibrated and validated using the SWAT-MODFLOW modeling. A total of 206 HRUs were generated for SWAT. A single-layer aquifer was simulated for the MODFLOW assuming that the connection between the shallow and deep wells is weak. Precipitation was set as the upper boundary for the SWAT-MODFLOW. Daily streamflow and monthly average groundwater level were used for the model's calibration and validation. Calibration and validation were performed using data obtained from 14 March 2013 to 31 December 2014 and from 1 January 2015 to 31 July 2015, respectively.

The input parameters of SWAT were processed and constrained within realistic uncertainty ranges based on the observed streamflow data. A sensitivity analysis was performed to identify the most sensitive parameters. In this study, 200 simulations were used for the sensitivity analysis. The sensitivities of the streamflow parameters are shown in Table 1. A small p-value indicates a lower probability of a nonzero regression coefficient. If the probability (p-value) is less than the significant value (0.05) of α, the null hypothesis is rejected, indicating that the results are statistically significant. Table 1 shows that the parameter most sensitive to streamflow is the runoff curve number, followed by the soil-saturated hydraulic conductivity and the soil evaporation compensation factor. The results indicated that geomorphology, soil properties, and temperature also affected the streamflow, but not significantly. MODFLOW was executed using the monthly data. Hydraulic conductivity was found to be the most sensitive parameter in groundwater modeling.

The water balance in the Wuling Watershed is presented in Table 2. Hydrological components, including the water yield, percolation, and evapotranspiration, were analyzed using 2-year data (1 July 2013–30 June 2015). The water yield represents the amount of water that leaves the watershed through reach, including the surface runoff, baseflow, and lateral flow (through the soil profile). The percentages of percolation, evapotranspiration, and water yield to precipitation were 10.94, 13.34, and 75.72%, respectively. Hence, three-fourths of the precipitation leaves the watershed through the reach, and approximately one-tenth of the total water percolates through the watershed. We observed that the year 2013 was wetter than usual and that the 2-year observation period (i.e., 2013–2015) was not long enough to comprehensively represent the weather in Wuling. Therefore, the 2-year observations were used to construct a 20-year baseline for the use of comparison with future projections.

The baseline data are required as a reference to simulate future scenarios (2021–2040). The 1986–2005 precipitation data for the Wuling Watershed were reconstructed based on the nearby long-term stations (i.e., the Hehuan, Ya-Kao, Tsui, Cui-Feng, and Zing-Kuan stations) using the Normal Ratio method. Meanwhile, the 1986–2005 temperature data for the Wuling Watershed were reconstructed based on the nearest meteorological station (i.e., the Sun Moon Lake Station) using an elevation-based correction method, that is, a decrease of 0.6 °C for every 100 m increase in elevation.

The RCPs form a set of greenhouse gas concentrations and emission pathways designed to support research on the impacts of and potential policy responses to climate change. The RCPs cover the range of the forcing levels associated with emission scenarios published in the literature. Four RCPs (i.e., 2.6, 4.5, 6.0, and 8.5) of the 21st-century pathways of greenhouse gas emissions, atmospheric concentrations, air pollution emissions, and land use were included in the Intergovernmental Panel on Climate Change (IPCC) fifth assessment report (AR5), as shown in Table 3. The downscaled precipitation and temperature projections under the RCP2.6, 4.5, 6.0, and 8.5 climate change scenarios were used to investigate the impact of climate change on major components of the hydrological cycle. Among the projections, RCP8.5 had the strongest driving forces for climate change.

Climate change induces variation in the water distribution, as shown in Table 4. The future precipitation projection was approximately 10.1% higher than the baseline precipitation projection. For the baseline water cycle components in the Wuling Watershed, percolation, evapotranspiration, and water yield accounted for 12.3, 15.2, and 72.5% of the total water, respectively. For the future water cycle components, percolation, evapotranspiration, and water yield accounted for 12.3, 14.3, and 73.4% of the total water, respectively. The percentage of percolation components remained the same, whereas evapotranspiration and water yield decreased and increased, respectively. Besides the changes in water distribution, the annual amounts of the studied hydrological components all increase. The percolation, evapotranspiration, and water yield increased by 10.9, 3.3, and 11.3%, respectively. This indicates that the hydrological cycle accelerated and the water distribution changed in the hydrological components.

Hydrological components are directly related to the available water resources. Our results showed that the hydrological components had different responses in the future projection as compared to the baseline projection. Precipitation increased in the wet seasons and decreased in the dry seasons under all RCP scenarios. The temperature has consistently increased throughout the year. Among the four RCP scenarios, RCP8.5 was the most drastic scenario, in which the changes in the percolation and streamflow corresponded highly to that of the precipitation. The lower streamflow in the dry season indicated less available surface water, while the higher streamflow in the wet season did not necessarily mean there was more surface water available. Instead, flood disasters and high-water turbidity could occur more frequently, and these conditions are not favorable within the context of water use. In particular, the amount of water increased in all water components and thus indicated an acceleration of the water cycle. In future climate projections, precipitation will tend toward white noise. Streamflow, evapotranspiration, and percolation became antipersistent, whereas the groundwater remained persistent. Thus, climate change significantly impacts surface water, while groundwater may serve as a buffer to mitigate the short-term impact on water resources, making groundwater fundamentally important for the water supply and other ecosystem-linked functions.

In this study, we proposed an approach that integrated physics-based meteorological and hydrological models to characterize and predict changes in hydrological components caused by climate change for a data-scarce mountainous watershed. Data limitations and model assumptions were embedded in the proposed approach. First, to ensure consistency with the groundwater level data resolution, hydrological data were converted from hourly to daily. Therefore, extreme hydrological events such as peak flow were not modeled at short intervals. Nevertheless, this did not prevent the necessity of using spectral analysis to assess the impact of climate change on 10-day, monthly, and seasonal trends by comparing the power law exponents of two equal-length periods for time series from the baseline and future scenarios. In addition, land use, topography, and soil distribution were assumed to remain the same during the SWAT modeling period. The IPCC sixth assessment report (AR6) can be used to further improve the modeling scenarios (IPCC 2022). Third, the slope effect on SWAT modeling is a concern (Bieger et al. 2015). Recent studies (Yacoub & Foguet 2013; Lu & Chiang 2019) have shown that SWAT is applicable for modeling mountainous watersheds (Yacoub & Foguet 2013; Lu & Chiang 2019). The denser HRUs help to delineate the hydrologic response to the topographic variation. In addition, the study site of the present study is small and has a mild topography, and the satisfactory calibration and verification results support the use of SWAT modeling.

The short-term hydrological data, physics-based meteorological and hydrological models, and spectrum analysis were integrated to evaluate the hydrological response to climate change for a data-scarce watershed. The Wuling Mountain Watershed in central Taiwan was taken as the study site. Twenty-eight-month hydrological data from 2013 to 2015 were used for the hydrological characterization. The results showed that the combined use of the physics-based hydrological and GCM models provided a feasible approach for characterizing the hydrological component distribution under the climate change stress for data-scarce watersheds. This integrated approach can be applied to other data-scarce watersheds to measure the impact of climate change on future water scarcity. The results of scenarios from CMIP5 for the climate projection for the period 2021–2040 showed that the future precipitation, water percolation, and streamflow will decrease by about 10% and increase by about 20–25% in the dry and wet seasons, respectively. The evapotranspiration is lower than that of the baseline in January, March, and December, whereas it can be as high as 4% during the other months. Under the RCP8.5 scenario, annual amounts of all studied water components increase. Percolation, evapotranspiration, and water yield change from 12.3, 15.2, and 72.5% of the total water for the baseline to 12.3, 14.3, and 73.4%, respectively. The increase of the water amounts and changes of the water distribution indicate the possible acceleration of the water cycle. Spectrum analysis showed that surface water tends toward becoming more antipersistent and harder to predict, while groundwater remains mildly persistent. Thus, groundwater may serve as a buffer for the impact of climate change on water resources for the Wuling Mountain Watershed.

This study was financially supported by the Ministry of Science and Technology of Taiwan under Grant No. MOST 111-2116-M-006-015 and MOST 111-2221-E-006-004. The data provided by the Central Geological Survey of Taiwan are greatly appreciated. The comments of the editor and two anonymous reviewers are greatly appreciated.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.