Assessing hydroclimate response to land use/cover change using coupled atmospheric-hydrological models

Abstract

Modeling techniques provide a straightforward means to dissect regional hydroclimate in response to changes in land use conditions. This study uses the Weather Research and Forecasting Model (WRF) and WRF-Hydrological modeling system (WRF-Hydro), driven by survey-based land use data in 1995 and 2015, to assess how central Taiwan’s hydroclimate responds to land use/cover change. We first run WRF-Hydro with observed rainfall as meteorological forcing to ensure reasonable runoff simulation, and then select ten cases under weak synoptic forcings in July and August in recent decades for the simulation under two land use conditions. The WRF-only simulation (i.e., uncoupled with WRF-Hydro) can reveal significant changes in heat fluxes and surface variables due to land use/cover change, including sensible and latent heat fluxes, 2-m temperature and specific humidity, and precipitation over the hotspots of urbanization or downwind areas. Coupling WRF with WRF-Hydro discloses varied runoff characteristics subject to land use/cover change: a general increase in average peak flow (~ 4.3%) and total runoff volume (~ 5.0%) accompanied by less definite time-to-peak flow, indicating a synergistic but sometimes competitive relationship between the pure hydrologic/hydraulic perspective and land–atmosphere interactions. By taking the difference between the uncoupled and coupled simulations, we verify that surface pressure, precipitation, and soil moisture are more sensitive to a better depiction of terrestrial hydrological processes; differences in the spatial variances of soil moisture can be as high as two orders of magnitude. Our findings highlight the importance of more comprehensive model physics in regional hydroclimate modeling.

Introduction

Even though land–atmosphere interactions involve complex physical processes and even human interventions (Pokhrel et al. 2017), scientific advancements have opened up an avenue for analyzing regional atmospheric circulations and water budget partitioning using atmospheric and/or hydrological models (Arnault et al. 2016; Kerandi et al. 2017). Through the use of various modeling techniques, numerous studies have pointed out that surface/subsurface conditions can exert considerable influence over the regional to global climate system (Seneviratne et al. 2010). For instance, Pal and Eltahir (2001) found that soil moisture (SM) plays a critical role in determining rainfall processes and the hydrologic cycle in the Midwest of the U.S. from late spring to early summer based on a regional climate model. A coupled subsurface, land surface, and atmospheric modeling study by Lo and Famiglietti (2011) indicated that the increase in groundwater, which leads to changes in SM content, evapotranspiration, surface temperature, and cloud cover, could be a determinant of uneven variations in global precipitation.

Among various drivers that modulate land surface conditions, land use/cover change (LUCC) that has been taking place over many regions of the world can explicitly influence land–atmosphere interactions through altering regional thermodynamic conditions (e.g., solar radiation, latent heat, and sensible heat), thereby yielding significant impact on local and regional climate (Salazar et al. 2015; Findell et al. 2017). In this regard, abundant modeling studies have also tried to disclose the relationship between hydroclimate and LUCC at varied spatial scales. LUCC examined in these studies can be classified into several forms, such as changes in vegetation cover (Cao et al. 2015; Qiu et al. 2021), irrigation patterns (Lo and Famiglietti 2013; Lo et al. 2021), deforestation/afforestation (Chen et al. 2019; Boysen et al. 2020), and urbanization (Chen et al. 2020; Hyka et al. 2022; Sharma et al. 2022). Regardless of the forms of LUCC, all these studies reached a similar conclusion that local and regional hydroclimate reflected through temperature, precipitation, or water cycle can respond to LUCC to a significant extent.

In Taiwan, there exists some valuable studies that employ atmospheric models to investigate the source of predictability (Fang et al. 2011), to develop a quantitative precipitation forecast system (Hong et al. 2015), as well as to explore how a hypothetical urban expansion influences local convection (Lin et al. 2008). In particular, Chen et al. (2020; C20 hereafter), as the predecessor of this work, applied the Weather Research and Forecasting model (WRF) to study how hydroclimate responds to LUCC in central Taiwan, where the most aggressive LUCC in the form of urbanization has been witnessed since the 1990s. One of the highlights of C20 was the use of realistic, survey-based land use data acquired from the National Land Surveying and Mapping Center (NLSC) in 1995 and 2007 to drive the WRF simulation. However, the NLSC released an updated version of land use data in 2015 that can better represent the present-day surface condition of Taiwan. Further, land–atmosphere interactions in C20 were primarily accounted for by the Noah land surface model (LSM, Livneh et al. 2011), which oversimplified some essential processes in terrestrial hydrology.

Although there have been numerous studies applying either atmospheric or hydrological models to investigate the impact of LUCC on regional hydroclimate, we still find a research opportunity to enrich the literature on using a coupled atmospheric-hydrological modeling framework for simultaneously discussing the responses of both atmosphere and land to LUCC, especially for an island with complex terrains such as Taiwan. Further, it would be essential to use the most realistic land surface data (i.e., land use and soil) for the modeling experiments. We thus set our study objectives as follows:

• 1. To conduct coupled atmospheric-hydrological modeling experiments to assess hydroclimate response to LUCC;

• 2. To assess the relative sensitivity of various hydroclimate variables to LUCC or different model physics.

In this study, as a successor to C20, we aim at (1) reinitiating a new batch of WRF simulations using the most updated land use data and (2) conducting another batch of simulations that couple WRF with the WRF Hydrological modeling system (WRF-Hydro). The first task will allow us to examine whether those previous findings of C20 remain valid in a scenario of continuous “urban sprawl,” whereas the second task will enable the discussion regarding variations in surface hydrologic processes (e.g., runoff) subject to LUCC. Moreover, we will be able to quantify how much the incorporation of more physics in hydrology into modeling can contribute to each field of hydroclimate variables by comparing coupled against uncoupled modeling results.

We briefly describe our assessment process next and append the flowchart (Additional file 1: Fig. S1) to this article. We commence with acquiring the necessary data for conducting our experiment, including land use, rainfall, and streamflow. To make the acquired land use data compatible with the modeling environment, we introduce some Geographic Information System (GIS) processing steps to the land use data. We will further describe the study area, data, and land use processing in "Materials". After selecting several cases under weak synoptic forcings for our simulations, WRF and WRF-Hydro must undergo some preprocessing prior to their runs. In addition, before we can couple WRF with WRF-Hydro to examine differences in simulated streamflow subject to LUCC, we perform a standalone run of WRF-Hydro (i.e., uncoupled with WRF) to examine whether it can produce reasonable streamflow with observed rainfall. We will elaborate more on the configurations of WRF and WRF-Hydro along with case selection in "Methodology". We are then able to conduct the WRF-only and coupled WRF-WRF-Hydro (WWH) simulations of the selected cases, driven by the processed land use in 1995 or 2015. We thus obtain four sets of simulation results for comparison, namely WRF with the land use in 1995 and 2015 (W95 and W15) and WWH with the land use in 1995 and 2015 (WWH95 and WWH15), for us to assess three types of modeling differences: (1) W15 – W95, as the baseline/uncoupled modeling response to LUCC; (2) WWH15 – WWH95, as the coupled modeling response to LUCC; and (3) WWH15 – W15, as the difference in model physics (by setting the land use as a control variable). We will present the three types of assessment results in "Results", and then discuss our findings in comparison with others in "Discussion". We will provide the concluding remarks and future research directions in "Conclusions and recommendations".

Materials

Study area

In consistency with C20, the study area of interest is central Taiwan, which includes five administrative regions (i.e., Miaoli, Taichung, Changhua, Yunlin, and Nantou counties/cities) with a total area of ~ 10,507 km². In this area, there are six major river basins including the Houlong, Dan-an, Dajia, Wu, Zhuoshui, and Beigang river basins from the North to South (Fig. 1). Most of these river basins originate from the Central Mountain Range, and all of them flow towards the West to the Taiwan Strait, part of the North Pacific Ocean.

Fig. 1

WRF’s three-level nested modeling domains and the area of interest (central Taiwan with the six major rivers)

Land use data

Not only did C20 indicate the significant variations in WRF simulations in response to different land use conditions but also the poor representation of Taiwan’s land use by the default land use data (e.g., the United States Geological Survey, USGS) in WRF (Cheng et al. 2013). To drive our newly designed WRF and WRF-Hydro simulations, we thus adopt Taiwan’s land use data released by the NLSC, which applied cadastral and field surveys, aerial images, and GIS analysis to develop the most realistic surface conditions in 1995 and 2015 (NLSC95 and NLSC15 hereafter). To make both the NLSC95 and NLSC15 data compatible with the WRF environment, we have to perform a series of GIS processing that reclassifies land use to match the ten major USGS land-use types (Fig. 2), and then converts/reprojects/resamples the data to the Lambert conformal projection at a resolution of 30 arc seconds (~ 1 km). This represents our model limitation. In Fig. 2, we can also note the dominant LUCC from the cropland and pasture type to the urban and built-up land, mostly over metropolitan and township areas in the five administrative regions. In the last step of land use data processing, we intentionally keep the NLSC95 data outside central Taiwan the same as the NLSC15 data to avoid our simulations being contaminated by LUCC in other areas of the island.

Fig. 2

NLSC land use data for central Taiwan in 1995 and 2015

Rainfall, streamflow, and meteorological forcing data

In addition to the land use data, we use rainfall, streamflow, and meteorological forcing data to either drive or assess the model simulations in this study. Regarding rainfall data, we acquire gauge and gridded rainfall data from the Central Weather Bureau (CWB) of Taiwan and the Taiwan Climate Change Projection Information and Adaptation Knowledge Platform (TCCIP), respectively. We use a total of 216 CWB gauges with hourly rainfall records for determining cases for WRF simulations. We use the TCCIP gridded data (1 km, daily), as the best available gauge-based precipitation data (Li et al. 2021) and the major meteorological forcing, for calibrating WRF-Hydro in a standalone mode. We will explain case selection and the calibration of WRF-Hydro more in the next section. Because the calibration of WRF-Hydro is based on the comparison between simulated and observed streamflow series, we acquire hourly streamflow data from Taiwan’s Water Resources Agency. The spatial distribution of stream gauges used is shown in Fig. 1. The standalone run of WRF-Hydro requires additional meteorological forcing data (i.e., incoming shortwave and longwave radiation, specific humidity, air temperature, surface pressure, and wind), which can be obtained from a preliminary coupled WWH simulation of each selected event for calibration.

Methodology

WRF configuration

We use WRF version 4.2.1 with the Noah LSM to conduct our simulations. As in C20, we run WRF with three nested domains (Fig. 1) at different horizontal resolutions, ranging from 18 km in the coarsest outermost domain (D1), 6 km in the middle (D2), and 2 km in the finest innermost domain (D3). There exists a stretched vertical grid consisting of 45 sigma levels. Regarding various parameterization schemes (e.g., cumulus, radiation, and cloud microphysics), we also opt for the same schemes as in C20. To supply WRF with necessary initial and boundary conditions, we use the National Centers for Environmental Prediction’s Final Reanalysis data. After the simulation of each selected case is completed, we take outputs from D3 for successive analysis. Please note that WRF based on our configuration has been thoroughly tested and verified with observations in our previous studies (e.g., C20 and Li et al. 2021); for instance, a sensitivity analysis has been performed to ensure a reliable configuration of WRF that utilizes the most adequate cumulous parameterization scheme and nested domains.

Apart from the basic configuration as described above, we modify the soil data in the WRF environment because Lin and Cheng (2016) have pointed out that the default soil texture in the WRF-Noah framework, developed by the Food and Agriculture Organization of the United Nations 16-category soil types, is not representative enough for territories outside the United States. We also believe that the most realistic geographical distribution of soil classifications is a key determinant of the accurate calculation of SM-atmosphere interactions. This is especially important for the coupled WWH simulation, which refines the delineation of hydrological processes at a spatial resolution much finer than WRF’s (see "WRF-Hydro configuration and pre-processing"). We thus follow Lin and Cheng (2016) to remap the local soil dataset, obtained from the Council of Agriculture, Executive Yuan, to the WRF soil classifications, and then replace the default WRF soil dataset with the local product. We also update the soil hydraulic parameters in the lookup table, including the slope of the retention curve, saturated matric potential, saturation hydraulic conductivity, porosity, field capacity, and wilting point. Please refer to Lin and Cheng (2016) for more details. Lastly, we should note that in assessing three types of modeling differences ("Introduction"), we keep the physical parameterization schemes and the soil data intact. In other words, the only thing that induces modeling differences is the change in land use. Such assessment can also eliminate systematic error in the numerical simulations due to an imperfect model (Bauer et al. 2015), thereby making WRF a valuable tool for unveiling the impact of LUCC on regional hydroclimate.

WRF-Hydro configuration and pre-processing

WRF-Hydro is a community distributed hydrological model able to simulate hydrologic responses, including surface runoff, subsurface flow, and channel flow. WRF-Hydro can be run in a standalone mode or coupled with WRF and an LSM (e.g., Noah) to provide feedback from a better delineation of surface/subsurface hydrological processes to the atmosphere (Gochis et al. 2020). The major physics behind WRF-Hydro entails several modeling components designed to enhance the calculation and spatial distribution of energy fluxes, SM, and overland and subsurface flows. Regarding flow routing, WRF-Hydro employs the diffusive wave equation and the classic Muskingum-Cunge approach for overland flow routing and channel flow routing, respectively. In this study, we use WRF-Hydro version 5.1.1 and adopt the HydroSHEDS elevation data (Lehner et al. 2008) at a resolution of 3 arc seconds (~ 90 m) for creating the routing grids in the pre-processing step. We determine the spatial resolution of WRF-Hydro routing grids as 100 m by setting the regridding factor as 20 (i.e., 2-km resolution in WRFs’ D3 divided by 20). We also determine the minimum contributing area as 2 km² by setting the number of routing grids to define the stream as 200.

Prior to running coupled WWH simulations, we calibrate WRF-Hydro and assess whether the hydrological model, based on the above configuration, can produce reasonable surface runoff subject to observed rainfall in a standalone mode. WRF-Hydro offers a wide range of parameters (up to several hundred) able to affect model outputs. For instance, there are three “table” files, namely SOILPARM.TBL, CHANPARM.TBL, and GENPARM.TBL, in which soil parameters (e.g., saturated water content and wilting point water content), Manning’s roughness coefficients, and soil hydraulic properties in the modeling domain can be found, respectively. Because of the modification of the soil data that have been made in the WRF environment as described in "WRF configuration", we believe the modified soil data, as well as the soil-related parameters, should be able to represent the most realistic spatial distribution in the modeling domain. We thus only perform the calibration of Manning coefficients (in CHANPARM.TBL) over different stream orders to achieve better correspondence between simulated and observed runoff series. Our manual calibration is based on a stepwise approach similar to Yucel et al. (2015). To be specific, we use a scaling factor to change Manning coefficient values within their valid physical range, and then determine the optimum parameter values by several performance measures (e.g., R² and root mean squared error, RMSE). As indicated in "Rainfall, streamflow, and meteorological forcing data", we adopt the TCCIP precipitation data as the major meteorological forcing in this calibration. Since hourly data are required for the meteorological forcing of WRF-Hydro, we perform temporal disaggregation of the daily TCCIP data by multiplying its daily rainfall amount with the rainfall percentage in each hour of a day, obtained from a CWB gauge identified with the greatest rainfall amount for a selected event. We obtain other necessary meteorological forcing data and antecedent SM states from the outputs of a preliminary coupled WWH run.

We adopt an event-based calibration/validation that includes two storm events on September 27, 2016, and July 30, 2017; the accumulated rainfall patterns corresponding to the two events are shown in Additional file 1: Fig. S2. We plot the calibration results in Additional file 1: Fig. S3, with several performance measures also denoted in the figure; the “satisfactory to very good” correspondence (Moriasi et al. 2015) between the simulated and observed runoff series suggests that the calibrated WRF-Hydro can produce reasonable hydrologic responses and is ready for the succeeding coupled modeling experiment. The calibration period of the hydrological model is expected to fall in line with the simulation period of each selected case, and the event-based calibration is common in copious studies relying on such computationally expensive hydrological models as WRF-Hydro (e.g., Yucel et al. 2015; Pasquier et al. 2022; Wang et al. 2022). If more computational resources were available, using automated calibration methods, including more parameters in the calibration (e.g., Silver et al. 2017), or performing a continuous and more extended simulation for the calibration could be another viable option.

Case selection and simulation time

In consistency with C20, we select cases under weak synoptic forcings for our modeling experiment to manifest the modeling differences resulting from LUCC. The idea is to circumvent substantial weather systems (e.g., typhoons or weather fronts) from outweighing land–atmosphere interactions. We hereby adopt the ten most representative convective rainfall events on a day in July and August from 1999 to 2016 (Additional file 1: Table S1), identified by C20, with a specific screening procedure that uses all the available CWB gauges in central Taiwan. Please refer to C20 for the details on the screening procedure.

For each selected case, the simulation begins at 0000 Coordinated Universal Time (UTC) one day ahead of the event (Day 0, 0800 local standard time, LST) and ends at 1800 UTC two days after (Day 3, 0200 LST), adding up a total simulation time of 66 h. The simulation time is longer than that in C20 because we have to account for the completeness of a hydrograph (e.g., from rising to recession limbs) that typically has a longer duration than the rainfall event. We take the first 12 simulation hours as the model spin-up period, and then obtain selected output variables from Day 1, 0000 to 2400 LST (except rainfall from Day 1, 1200 to 1800 LST, for analyzing convective rainfall only in the afternoon) for the follow-up differential analysis. Regarding surface runoff, we obtain the output information from Day 0, 2100 to Day 2, 2300 LST, adding up a total runoff duration of 50 h.

Results

WRF simulation under LUCC (W15 – W95)

To assess how WRF responds to different land use conditions, we take the difference between the mean field of the ten selected cases driven by NLSC15 and that driven by NLSC95. Variables of interest include sensible and latent heat fluxes (SH and LH), 2-m temperature and specific humidity (Q), planetary boundary layer (PBL) height, surface pressure, 10-m wind, and precipitation, and the results pertaining to W15 – W95 are shown in Additional file 1: Figs. S4 and S5. The figures show a clear increase (decrease) in SH and temperature (LH and Q) over the western part of the study area, where urban expansion is evident. In addition to the urban heat island effect, we believe that the notable reduction in irrigation activity and surface/subsurface runoff, as a result of the transformation from the irrigated cropland into the urban and built-up land (Fig. 2), should be another determinant of the changes in heat fluxes and temperature. The PBL pattern bears a high similarity with the temperature since the increase in the thickness of PBL results from more enhanced atmospheric mixing by a warmer surface. The diminished surface pressure is also related to the warmer surface, which causes a drop in air density, favoring the upward motion of air. In association with the PBL and pressure patterns, local convergence of 10-m winds and intensified onshore vectors can be identified. At last, precipitation that symbolizes an integrated response to the above variables as well as local topography exhibits an intensification over the windward side of the hills rather than the hotspots of urbanization areas. More discussions regarding how LUCC can induce off-site influence on precipitation can be found in C20, which showed a series of tests for statistical significance (Benjamini and Hochberg 1995) and an explanation of physical mechanisms.

The above findings are basically in line with C20, so we can reasonably attribute the similar variations in the simulated fields to the resemblance between the land use in 2007 and that in 2015 used in C20 and the present study, respectively. In fact, a preliminary examination of the land use data in different years by C20 has already pointed out that the change in urban and built-up lands from 2007 to 2015 was not as dramatic as the change from 1995 to 2007, suggesting that any variations in the simulated fields driven by urban sprawl should have developed since 2007. C20 also indicated an increase in irrigated cropland areas (quite notable in Changhua and Yunlin) from 2007 to 2015, and we have corroborated that such an increase can play a role in mitigating the magnitude of the variations in local hydroclimate in response to urbanization in the present study. This finding implies that more optimized land management strategies and geoengineering could create a reversible scenario for adverse hydroclimatic conditions (Mostamandi et al. 2022).

Coupled WRF-WRF-Hydro simulation under LUCC (WWH15 – WWH95)

According to the set objective, the second task is to assess the coupled modeling response to LUCC; that is, we derive the mean field of the ten WWH simulations driven by either NLSC15 or NLSC95, and then calculate the difference between the two mean fields, referred to as WWH15 – WWH95. As shown in Additional file 1: Figs. S6 and S7, the results of coupled modeling bear a high resemblance to the results of uncoupled modeling (Additional file 1: Figs. S4 and S5) in general. The spatial patterns and magnitude of heat fluxes, temperature, and PBL are very similar between the two modeling experiments, suggesting that the variations in these variables are determined more by LUCC and less sensitive to the add-on configuration of terrestrial hydrological processes and physics. In other words, an enhanced representation of water bodies and river routing seems to yield less influence over the above variables. In contrast to the above variables, surface pressure and precipitation exhibit more distinct features in the magnitude or spatial patterns between the two modeling experiments. We will further examine and discuss whether the varied sensitivity of each variable to the more comprehensive model physics makes sense in the following subsection.

One of the merits of running the coupled WWH model is being able to analyze runoff characteristics subject to LUCC. In each of the ten cases, we identify or derive three runoff statistics, namely average peak flow (Q_p), time-to-peak flow (T_p), and total runoff volume (V) in the 50-h duration, over the six watershed outlets (green dots in Fig. 1) for runoff analysis. To exclude the effect of spatial scale from our analysis, we divide Q_p and V by the drainage area of each watershed before taking the watershed-wide average. We are then able to calculate the difference between each runoff statistic in WWH95 and that in WWH15, as shown in Table 1. The overall average of the differences over the ten cases discloses a universal increase in every runoff statistic under LUCC: Q_p, T_p, and V all increase by approximately 4.3%, 0.5 h, and 5.0%, respectively. Even though we find the universal increase in the runoff statistics on average, we notice that the differences in T_p exhibit a prominent inter-case disagreement, with half of the cases being greater or less than 0 h. We will further discuss the impact of LUCC on runoff statistics in "Impact of LUCC on runoff statistics based on coupled modeling".

Table 1 Runoff statistics in central Taiwan derived from each WWH95 or WWH15 case simulation

Uncoupled vs. coupled simulation under the same land use (WWH15 – W15)

To further examine which variables are more sensitive to the model physics, we use and control the land use in 2015, and then calculate the difference between the mean field in the WWH15 simulation and that in the W15 simulation (i.e., WWH15 – W15). Figure 3 shows that, in comparison with the previous results (Additional file 1: Figs. S4 or S6), SH and temperature are less sensitive to the add-on configuration of terrestrial hydrological processes and physics. On the other hand, we can find a more evident increase in LH and Q, suggesting a wetter overland scenario in the coupled WWH simulation. While previously claiming that the difference in PBL in the uncoupled simulation resembled much to that in the coupled simulation (Additional file 1: Figs. S5 vs. S7), we find a peculiar pattern of PBL in Fig. 4a, as a result of the subtle difference in atmospheric mixing induced by the slight contrast in temperature (Fig. 3c). Surface pressure, 10-m wind, and precipitation are three variables showing the most significant modeling differences in the uncoupled vs. coupled simulation. The differences in these three variables are actually comparable to or even larger than the LUCC-induced differences, which can be corroborated by Table 2, where we show the spatial variance of each simulated variable in central Taiwan derived from three types of modeling differences.

Fig. 3

Difference between the mean field (i.e., average over ten selected cases) in the WRF-WRF-Hydro simulation and that in the WRF-only simulation (both driven by the NLSC15 land use, i.e., WWH15 – W15), with respect to (a) sensible heat, (b) latent heat, (c) 2-m temperature, and (d) 2-m specific humidity

Fig. 4

As in Fig. 3, but with respect to (a) planetary boundary layer height, (b) surface pressure, (c) 10-m wind, and (d) precipitation

Table 2 Variance of each simulated variable in central Taiwan derived from three types of modeling differences

The above results demonstrate that each simulated variable displays a varying degree of responses to the change in land use or model configuration, and some of the responses are highly associated with the extent of the moist condition near the terrestrial surface. We thus take one step further to examine the top-layer SM from each simulation set by plotting its differences according to three types of modeling differences in Fig. 5. In the first two modeling differences, we can see a slight increase (decrease) in SM over the cropland areas adjacent to the cities/townships (mountainous regions). In contrast, the third modeling difference shows that the magnitude of change in SM due to the use of WRF-Hydro is much greater than that due to LUCC; the increase in SM is immense from upstream to downstream areas except in some mountainous regions. In fact, the spatial variances of SM verify that the difference can be up to almost two orders of magnitude (Table 2). We will further rationalize the increase in SM and discuss the relative sensitivity of atmosphere and land variables to LUCC or model physics in "Relative sensitivity of atmosphere and land variables to different forcings".

Fig. 5

Difference in the mean top-layer soil moisture according to three types of modeling differences assessed in this study

Discussion

Impact of LUCC on runoff statistics based on coupled modeling

The result of runoff analysis suggests that the intensification of rainfall over the windward side of the hills (i.e., upstream areas) due to LUCC can generate higher peak flow and more runoff on the surface. Still, urban sprawl that takes place over the mid-to-downstream of the watersheds yields indefinite influence on the time when the peak flow arrives at the downstream outlet. This result is intriguing because the supposedly increasing flow velocity (i.e., shortened T_p) due to the expansion of built-up or “impervious” areas is not always valid. The indefinite influence on T_p also implies that the pure hydrologic/hydraulic point of view (i.e., lower surface roughness, thus higher flow velocity) might not be identified in reality owing to the existence of its offset against intricate land–atmosphere interactions.

Even though many studies have made similar attempts to assess the impact of LUCC on runoff statistics, most of these studies drew their conclusions based on using solely hydrological models (e.g., Hu et al. 2021; Liu et al. 2023), uncoupled simulations (e.g., Pasquier et al. 2022), or coupled earth system models with land surface models at a coarser spatial resolution (e.g., Qiu et al. 2021). Our finding pertaining to the competitive relationship between the pure hydrologic/hydraulic perspective and land–atmosphere interactions is thus unique for using the coupled atmospheric-hydrological model at an island scale and a finer resolution. Nevertheless, it should be noted that our finding is subject to some of our method limitations, such as the processed land use data at a resolution of 30 arc seconds and the selected cases under weak synoptic forcings.

Relative sensitivity of atmosphere and land variables to different forcings

The pattern of increased SM, derived from the difference between the uncoupled and coupled simulations, is deemed legitimate since the higher SM mostly takes place along or nearby the more realistic delineation of stream links (Fig. 1). This finding indicates that SM is most sensitive to the change in model physics. More frequent infiltration and higher SM were also found by Wang et al. (2020) in their WWH simulation compared against the WRF-only simulation of storm events. The higher SM not only explains the wetter overland scenario (thereby promoting the higher Q and LH) but also illustrates the importance of a more physical representation of water bodies and river routing at a finer resolution in atmospheric modeling for addressing simulation differences subject to various forcings.

For each simulated variable in Table 2, we highlight the highest variance among the three modeling differences to reveal the relative sensitivity of each variable to the three types of modeling differences. We find that heat fluxes, temperature, Q, and PBL are more sensitive to the first or second type of modeling difference, suggesting LUCC can exert more influence on regional hydroclimate than the change in model physics. In contrast, surface pressure, precipitation, and SM show higher sensitivity to the change in model physics than LUCC. The relative sensitivity of each simulated variable, reflected through the comparison among spatial variances, also implies the level of accuracy for a specific variable that could be attained in the context of hydroclimate prediction. For instance, coupling WRF with WRF-Hydro can induce a ~ 30% increase in the spatial variance of the precipitation field, which can dramatically impact the accuracy of precipitation prediction regarding both spatial patterns and amounts. This finding again illustrates the advantage and necessity of accounting for more comprehensive model physics in regional hydroclimate modeling (Naabil et al. 2023).

Conclusions and recommendations

This study obtained Taiwan’s survey-based land use data in 1995 and 2015 to assess central Taiwan’s hydroclimate response to LUCC. The assessment was conducted by simulating ten selected cases under weak synoptic forcings in July and August using WRF. To enhance the depiction of terrestrial hydrological processes, WRF was further coupled with WRF-Hydro; how the coupled simulation differed from the uncoupled simulation was then examined. Our major findings and concluding remarks are summarized as follows:

• 1. In response to urbanization, as the most notable LUCC of central Taiwan in the last two decades, the WRF-only simulation can generate reasonable modeling differences, including increased (decreased) H and temperature (LE and Q) over the hotspots of urbanization areas. The elevated PBL, along with enhanced local convergence and onshore wind, can induce intensified rainfall over the downwind areas. The adverse hydroclimate response (e.g., urban heat island) became milder in comparison with C20, partly associated with the increased irrigated cropland areas from 2007 to 2015.

• 2. The results of coupled modeling (i.e., WWH) generally bear a high resemblance to the results of uncoupled modeling (i.e., WRF-only) in response to LUCC, yet we found some variables (e.g., SM) being highly sensitive to the add-on configuration of terrestrial hydrological processes and physics. Moreover, coupled modeling enabled the analysis of runoff characteristics, of which we found a general increase in average peak flow and total runoff volume in correspondence with the upstream intensification of rainfall due to LUCC. On the other hand, time-to-peak flow exhibited a prominent inter-case disagreement, with half of the cases being delayed or expedited.

• 3. The less definite influence of LUCC on time-to-peak flow implies a synergistic but sometimes competitive relationship between the pure hydrologic/hydraulic perspective and land–atmosphere interactions. This unique finding highlights the importance of using the coupled atmospheric-hydrological model for an island-scale, fine-resolution modeling study.

• 4. The difference between the uncoupled and coupled simulations corroborated the more sensitive response of certain variables (e.g., surface pressure, precipitation, and SM) to the better depiction of terrestrial hydrological processes than LUCC. In fact, WWH can provide an explained variance to a much more significant extent than WRF-only for SM, illustrating the necessity of accounting for more comprehensive model physics in regional hydroclimate modeling.

Since using the most updated land use data plays a vital role in such LUCC studies as the present one, we are in the process of obtaining the most updated version of the survey-based land use data from the NLSC to facilitate our future research for Taiwan. As long as the influence of substantial weather systems on regional hydroclimate can be separated from that of LUCC, another future research opportunity is to adopt a continuous and more extended simulation compared to the case-based assessment in the present study. Last but not least, assessing the combinatory impact of LUCC, model physics, and soil texture on hydroclimate modeling also merits further investigation.