Seasonality and Impact Factor Analysis of Streamflow Sensitivity to Climate Change Across China

Streamflow sensitivity to climate change is an important indicator for evaluating the effects of climate change on terrestrial water. This study analyzed the spatial pattern and seasonality of streamflow sensitivity to climate change for 425 catchments across China. The results indicate that precipitation is consistently a more important contributor than temperature to streamflow variability (mean streamflow sensitivity to annual precipitation [εP,a] = 0.92; absolute mean sensitivity to annual temperature [εT,a] = 0.05). Meanwhile, the seasonal sensitivity evaluation found that streamflow response to precipitation change in warm seasons (mean εP,w = 0.90) is significantly greater than in cool seasons (mean εP,c = 0.46), but the sensitivity of streamflow to changes in temperature was stronger in cool seasons than in warm seasons. The magnitude of streamflow sensitivity to temperature in cool seasons is significantly affected by elevation (Pearson's r = 0.31, P < 0.01), which is possibly due to the increased snow and glacial melt volume of the high‐elevation catchments in cool seasons. We also found that catchments with higher elevation, drier climate, less vegetation cover, and lower infiltration rates are more likely to present seasonally asymmetric precipitation sensitivity. A possible reason for this spatial pattern is that the highland and continental climate regimes generally exhibit strong precipitation seasonality, thus enhancing the asymmetry of streamflow response to seasonal precipitation change. These findings are essential to improve understanding of the impact of climate change on hydrological processes and thus for formulating adaptive water management strategies for future climate change in China.

. For example, Labat et al. (2004) reported a 4% global runoff increase in response to a 1°C global temperature rise based on a statistical wavelet-based method. For much finer information, analysis of streamflow sensitivity to climate changes in both hydroclimate annual means and seasonal variations can better inform assessments of societal and ecological vulnerability for future water resources (Konapala et al., 2020). For example, warming in cool seasons causes slower snow accumulation and earlier snowmelt, while warming in warm seasons may cause larger peak flows, which will strain reservoir capacities, and larger evaporative losses in summer (Li et al., 2017). Previous studies have adequately assessed streamflow sensitivity to climate change in terms of mean annual precipitation (Sankarasubramanian et al., 2001;Tang et al., 2019;Wang, 2020), temperature (Fu et al., 2007;, and other variables (Berghuijs et al., 2017;Xiao et al., 2020). Only a few studies have examined the distinct impacts of seasonal climate change on streamflow sensitivity, for example, in the Colorado River basin (Vano & Lettenmaier, 2014), the North American Pacific Northwest (Vano et al., 2015), Alaska (Curran & Biles, 2021), and the Western United States (Ban et al., 2020). The spatiotemporal distribution and seasonal patterns of streamflow sensitivity to climate change are largely unknown, particularly for China, where studies only provide information at a limited regional basin scale (Wang, 2020;Zhang et al., 2018) and thus provide little insight into diverse national climate sensitivity patterns through different hydroclimatic zones. Evaluating the streamflow sensitivity to climate change, especially at seasonal scales, is considered important to support water resources planning and management and to reduce the risks of water scarcity and water-related hazards in different seasons, because this information could lead to more seasonally specific recommendations on how best to adapt to future changes.
A variety of methods, including hydrological modeling (Ban & Lettenmaier, 2022;Vano et al., 2015), multivariate statistical regression (Vogel et al., 1999;Zhang et al., 2018), and the Budyko framework (Wang, 2020;Zheng et al., 2021), have long been under development to study streamflow sensitivity in changing environments. Hydrological modeling has generally been applied to simulate streamflow response to different forcing under specific climate change scenarios (Ban et al., 2020). But modeling approaches are somewhat unsatisfactory due to the impossibility of describing scenarios for the full, wide range of complex climate variables . For example, many studies use simplified temperature or precipitation addition approaches to construct warming or wetting scenarios (Das et al., 2011;Marshall et al., 2019;Singh & Bengtsson, 2004), which might ignore the internal correlations between climate variables and thus result in unreliable sensitivity results. Furthermore, climate sensitivity analyses performed on the same basin using different conceptual watershed models can lead to significantly different results (Sankarasubramanian et al., 2001). Multivariate statistical regression and Budyko framework are statistically robust and theoretically reliable and are widely applied to estimate the relationship between climate and streamflow; but they rely heavily on historically observed streamflow and climate change information for the basins of interest. Therefore, the reliability of streamflow and climate data sets plays a role in the climate sensitivity assessment of hydrological processes, and it also plays a key role in calibrating and validating climate and hydrological models Switanek et al., 2017).
With the rapid economic development in China, observed flow regimes have been greatly affected by human activities, including water diversion, water withdrawal, and reservoir operation (Gou et al., 2020;Miao et al., 2022). For instance, hydrological processes are highly regulated (flood peaks attenuated and low water increased) by upstream reservoirs after water storage (Xu, 1996). For example, China had built about 98,822 reservoirs by 2018, with a total storage capacity of about 895.3 km 3 (Tian et al., 2021). Relying on observed streamflow that is affected by human interference might mislead us as to the interpretation of streamflow sensitivity to climate change because such data contain major flow alterations due to human diversions and water management (Lehner et al., 2019;Miao et al., 2022). Distributed hydrological models that represent the natural rainfall-runoff relationship and the upstream-downstream flow propagation process while excluding the effects of human activities are a commonly used naturalization tool in hydrological studies (Terrier et al., 2021). Many past studies have used model-reconstructed natural streamflow data to assess climate change-related hydrological responses (Dariane & Pouryafar, 2021;Fu et al., 2007;Wang, 2020).
China's climate, which has greatly warmed over the past few decades (Piao et al., 2010;Xu, 1996), is extremely vulnerable to global warming due to influences from the East Asian monsoon and complex topography (Zhang & Zhou, 2020). Continued warming may lead to a further acceleration of the global water cycle, increase the frequency and intensity of hydrological extreme events such as floods and droughts (Jones et al., 2012), and even increase climate warming through water vapor feedback (Huntington, 2006). In addition, previous studies have 10.1029/2022EF003062 3 of 17 consistently projected increased surface mean temperature (You et al., 2021), increased evapotranspiration (You et al., 2021), and more frequent extreme climate events (Zhu et al., 2021) across China under Coupled Model Intercomparison Project Phase 6 simulations. An earlier study also found that warming in China is not evenly distributed on a sub-annual scale, with four times greater warming in the winter season than in the summer season across most of the region (Piao et al., 2010). Therefore, it is important to explore the regional and seasonal patterns of streamflow sensitivity to climate change, as well as to reveal the main impact factors that may affect adaptation to future climate change.
Given the above, the questions we address here are as follows: (a) How does the sensitivity of streamflow to changing climate differ across China? (b) How does the seasonal pattern of the streamflow sensitivity to changing climate vary? And (c) What are the main impact factors affecting the response of streamflow to annual and seasonal climate change? We address these questions by first providing a technical reconstruction of long-term natural streamflow for 425 catchments across 10 major river basins of China using a hydrological model and a post-processing bias-correction algorithm (Sections 2.2 and 2.3). We then discuss our development of a multiple regression statistical model to explore the relative sensitivity of streamflow to changes in annual and seasonal climate (Section 2.4). Finally, to assess the possible impact factors affecting streamflow sensitivity to climate change, we analyze the relationships between sensitivity and catchment characteristics (Section 2.4).

Data Sets
In this study, four kinds of data sets were used to explore seasonality and the possible impact affecting seasonal streamflow sensitivity to climate change across China: (a) forcing data sets for driving hydrological models; (b) natural (or near-natural) streamflow data sets for post-processing bias correction; (c) climate data sets for sensitivity analysis; and (d) catchment characteristics data sets for the impact factor analysis of streamflow sensitivity to climate change. Model forcing data sets include lateral inflows and river hydrography across China. The China Natural Runoff Data set version 1.0 (CNRD v1.0), which uses previously constructed 0.25° × 0.25° daily runoff data sets from a calibrated land surface model, was used to provide lateral natural inflows in this study. Details of the CNRD v1.0 quality-control procedures and model input data sets (e.g., vegetation, soil, and topographical data) are given in our previous work . A global vector-based river network data set, the MERIT-Basin (Multi-Error-Removed-Improved-Terrain-Basin) hydrography data set (Lin et al., 2019), was used to produce the flow paths for sub-catchment-scale streamflow simulation. MERIT-Basin provided sub-catchment domain and river centerlines for each river reach, along with the corresponding river reach properties, including catchment area, Strahler stream order, channel length, and upstream/downstream links.
Streamflow data sets of 425 hydrological stations ( Figure 1) were collected from the hydrological yearbook of China and local water resources departments. Among the 425 sets of station records, 189 gauge records are naturalized by the Ministry of Water Resources of China based on the water balance principle. The remaining stations selected have no upstream dams present or have a low level of dam influence. Thus the streamflow data is considered to be near-natural, that is, influenced relatively little by anthropogenic activity (e.g., dams and irrigation) with respect to the water balance, and suitable for long-term natural streamflow reconstruction. For simplicity, we regard these streamflow series from the 425 stations as "natural streamflow" hereafter. Natural streamflow records for the 425 gauge stations are available at a monthly time step and cover at least 10 years during the period 1961-1979 (this period chosen to avoid human activity interference that occured in China during its boom in socioeconomic activity after 1980). Drainage areas for the 425 gauged catchments are from 17.5 to 1,705,383 km 2 ( Figure S1 in Supporting Information S1), and span a range of hydroclimatic regimes, from arid (e.g., Northwest River drainage system), to semi-arid (e.g., Yellow River), to moderate (e.g., Huai River), to wet (e.g., Pearl River).
Climate data sets, including daily temperature, daily maximum and minimum temperature, and precipitation variables were obtained from the China Meteorological Administration (CMA) and the China Meteorological Forcing Data set (CMFD), which together covered the full research period from 1961 to 2018. CMA climate data are gridded and interpolated from ∼2,400 weather stations and cover most of the years 1961-2014 (Miao et al., 2022). The CMFD is a daily gridded meteorological forcing data set, available at 0.1° × 0.1° grid resolution during the period from 1979 to 2018, which comprises a fusion of remote sensing products, reanalysis data sets, and in situ station data (He et al., 2020). We carried out cross-validation tests for an overlap period  of the two meteorological data sets in previous work, and the results showed that the two data sets passed the consistency test . Potential evapotranspiration was calculated by the method of Hargreaves and Samani (1985) based on the daily minimum, maximum, and mean temperatures and extra-terrestrial radiation (Text 1 in Supporting Information S1).
We collected catchment characteristics for each drainage basin, including elevation, vegetation cover, field capacity, and saturated hydraulic conductivity. The 1 km digital elevation model (DEM) data are obtained from the Cold and Arid Regions Science Data Center at Lanzhou. The vegetation cover data were obtained from the Advanced Very High Resolution Radiometer (AVHRR)-based Global Inventory Modeling and Mapping Studies normalized difference vegetation index (NDVI) data set. We averaged the annual maximum NDVI values (range from 0 to 1) during the period from 1982 to 2015 to represent the general vegetation condition for each of the 425 catchments. The soil texture properties including field capacity and saturated hydraulic conductivity were obtained from the 30 arc s (∼1 km) China Data set of Soil Properties by Shangguan et al. (2013) and the China Data set of Soil Hydraulic Parameters by Dai et al. (2013), respectively.

Hydrological Model-Based Natural Streamflow Reconstruction
Compared to the observed streamflow, the reconstructed natural streamflow is more interpretable with respect to climate change sensitivity because it removes major flow alterations due to human diversions and water management. The Routing Application for Parallel computation of Discharge (RAPID) model version 3.13.0 (David et al., 2011) was used to reconstruct long-term  streamflow under natural conditions for 425 selected catchments. The RAPID model was designed with a vector river reach structure, which has a high level of efficiency and flexibility compared to traditional grid-based flow-routing models, such as the Lohmann routing model (Lohmann et al., 1996) and the Soil and Water Assessment Tool model (Arnold et al., 1998). The RAPID model supports matrix-based forms of the Muskingum method and is solved using a highly efficient linear system solver for reach-by-reach channel routing. Following Lin et al. (2019), the flood wave travel parameter k and the dimensionless weighting factor x were estimated for each vector river flowline. Lateral inflow Figure 1. Locations of the 425 hydrological stations used in this study. Colored relief shows relative topographic highlands (blue shading) and lowlands (white shading) using the 1 km digital elevation model data set. Gray boundaries indicate the 10 major river basin areas over China: I-Songhua River; II-Liao River; III-Hai River; IV-Yellow River; V-Huai River; VI-Yangtze River; VII-Southeast River drainage system; VIII-Pearl River; IX-Southwest River drainage system; and X-Northwest River drainage system. was produced by means of an area-weighted coupling interface supported by the ArcGIS-based Python toolbox (https://github.com/Esri/python-toolbox-for-rapid). We ran the RAPID model at a daily time step. In this study, the MERIT-Basin data set from Lin et al. (2019) and the System for Automated Geoscientific Analyses (SAGA) tool version 8.2.0 (Conrad et al., 2015) were used to automatically outline and aggregate the catchment boundaries from upstream rivers for the 425 catchments across China.

Statistics-Based Post-Processing Bias Correction
Raw streamflow simulations inevitably suffer from having a series of uncertainty sources because the river discharge model sits at the end of a hydrologic modeling chain consisting of atmospheric forcing, land surface models, and a river routing model (Lin et al., 2019). With the growing demand for the representation of natural hydrological processes and for exploring their responses to climate, higher-quality reconstructed streamflow data are being required (Miao et al., 2022). A statistical post-processing procedure can be used to remove the systematic bias and achieve reliable long-term natural streamflow reconstructions for hydroclimatology studies (Li et al., 2022).
The scaled distribution mapping (SDM) bias-correction method (Switanek et al., 2017) was used in this study to reduce streamflow simulation uncertainty while preserving raw hydrological change signals. The SDM method scales the referenced streamflow distribution (natural streamflow records at 425 stations during the period of 1961-1979 in this study) in magnitude by modeling simulated streamflow changes and the likelihood of events under nonstationary conditions. The scaling changes as a function of the bias-correction period, that is, 1961-2018. In this study, SDM scales the distribution of streamflow by a multiplicative or relative amount, and only values of positive streamflow exceeding a specified threshold (0.1 m 3 /s) are used to build the distributions. The SDM procedure is as follows. First, we use a gamma distribution to fit the baseline referenced natural streamflow (ref_data), baseline simulated streamflow (model_data), and long-term simulated streamflow (sce_data) values by the maximum likelihood method. The gamma distribution is expressed as where γ (>0) is a shape parameter, θ (>0) is a scale parameter, a (>0) is the variable of interest (streamflow in this study), and Γ(γ) is the gamma function evaluated at γ.
Next, we use the fitted γ and θ parameters to find the corresponding cumulative distribution function (CDF) values of the variable in the foregoing streamflow time series: ref_data, model_data, and sce_data. We then calculate the scaling between the fitted raw model_data distribution and the fitted raw sce_data distribution at all CDF values corresponding to the long-term sce_data time series. The scaling is calculated as where SF R is an array of relative scaling factors, ICDF sce and ICDF mod are the inverse cumulative distribution functions (ICDFs) for the fitted sce_data and raw model_data distribution, and CDF sce values are the estimated CDF values for the sce_data corresponding to the fitted distribution. We calculate recurrence intervals (RIs) for the three sorted streamflow arrays. The recurrence interval array is where CDF is the array of values found by Equation 2. The scaled (or adjusted) RI for the raw model_data is determined from where RI sce is the RI for the raw sce_data values and RI ref and RI model are linearly interpolated RIs for the ref_data and raw model_data, respectively. The maximum value, which is greater than or equal to 1, is used to evaluate each value in the RI scaled array. The corresponding scaled CDF values are obtained from where the CDF scaled array is subsequently sorted in descending order and reflects the scaling of the model_data change in event likelihood with respect to the ref_data likelihoods. The initial scenario array of bias-corrected values is calculated as where ICDF ref is the inverse cumulative distribution function for the ref_data fitted distribution, and CDF scaled and SF R are obtained from Equations 2 and 5.
For each catchment, we automatically corrected the simulated streamflow using the scaled fitted gamma distribution for CDF values corresponding to the streamflow time series. Four skill scores-Pearson's correlation coefficient (CC), percent bias (Pbias, unit: %), Nash-Sutcliffe efficiency coefficient (NSE), and Kling-Gupta efficiency coefficient (KGE)-were selected to explore the performance of natural streamflow reconstruction. The detailed calculation procedures for those skill scores can be seen in our previous work (Miao et al., 2022).

Streamflow Sensitivity to Climate Change
To explore streamflow sensitivity to climate change, the long-term catchment average monthly values of the climatic and hydrologic variables were calculated and analyzed. A simple linear statistical model is considered because it can approximate climate and streamflow data, and theoretical expressions are easily obtained for the streamflow sensitivity to changes in temperature ε T (T, Q) and changes in precipitation ε P (P, Q). A multivariate linear model relating Q to P and T is where ε T,w and ε P,w are model regression coefficients, and the second subscript for this coefficient, for example, ε T,w and ε P,w , indicates the period of calculation (a: annual, w: warm season, and c: cool season); Q, T, and P are the streamflow, temperature, and precipitation values for the specific calculation period; and ω is a residual. Note that T, P, and Q are all normalized to represent the relative streamflow sensitivity to climate change. For the annual research period, we used water year in all cases, where a water year is defined as the period from 1 May to 30 April of the following year and named according to the calendar year in which the period ends. Warm and cool seasons are defined as the periods from 1 May to 31 October and from 1 November to 30 April of the following year.
We used the bivariate ordinary least squares regression method to estimate the ε T,a , ε P,a , ε T,w , ε P,w , ε T,c , and ε P,c coefficients in Equation 7 based on statistical inputs for different periods, that is, a: annual, w: warm season, and c: cool season. For regression models, we conducted a Student's t-test to assess whether the streamflow was sensitive to climate variable at the 10% significance level. As in Ban et al. (2020), we used four indexes (Tperf, Pperf, Tsub, and Psub) to further quantify the relative strength of streamflow sensitivity to climate change during warm versus cool seasons. Taking temperature as an example, the calculation of Tperf and Tsub is as follows: Here, ε T,w and ε T,c are the streamflow sensitivity to temperature in warm and cool seasons, respectively. Tperf > 0 indicates that streamflow sensitivity to temperature change in a basin has the same response direction, that is, all positive or negative in warm/cool seasons, and vice versa. Tsub > 0 indicates that a basin has a larger streamflow sensitivity to temperature change in warm seasons than cool seasons and vice versa. Pperf and Psub use definitions analogous to those of Tperf and Tsub but for precipitation change.
To explore the potential influence of catchment characteristics on streamflow sensitivity to climate change, we analyze the correlation relationships between six characteristics and ε T,a , ε P,a , Tsub, and Psub indexes. The six catchment characteristics include aridity index (P/PET), runoff efficiency (Q/P), elevation (DEM, unit: m), vegetation cover (NDVI), field capacity (W, unit: mm), and saturated hydraulic conductivity (Ks, unit: mm/day).

of 17
Average catchment values of the six characteristics were calculated based on the database in Section 2.1 and the catchment boundaries generated in Section 2.2. Figure 2 shows natural streamflow reconstruction performance for the 425 hydrological stations, which is estimated using hydrological models with post-processing bias correction (Sections 2.2 and 2.3). Figure 2a shows the spatial distribution of reconstructed mean annual streamflow during the period 1961-2018 across the 425 catchments. The average streamflow for a water year ranges from near zero (e.g., in the arid region) to over 1,000 mm. Overall, the reconstruction performed well in representing the spatial pattern of streamflow as it decreases from the southeastern coast to northwestern lands in China. Similarly, basin-aggregated precipitation is generally high in humid southeastern China and low in northwestern China ( Figure S2a in Supporting Information S1). Figure 2b shows the relationship between the model area and the compiled observed drainage area (collected from the hydrological yearbook, local water resources department, and previous literature) for the 425 catchments. The total aggregated area over the entire 425-catchment domain is 45.75 × 10 5 km 2 (mean: 0.6 × 10 5 km 2 , standard deviation: 1.74 × 10 5 km 2 ). Note that the total area does not include nested catchments.

Natural Streamflow Reconstruction and Climate Change Analysis
In Figure 2b, we found that the modeled drainage area is closely matched to the collected observed drainage areas (Pearson's r = 0.99, P < 0.01, and mean absolute Pbias less than 10%), indicating the river networks of basin catchments are satisfactory within the RAPID river-routing model.  Figure 2d illustrates the probability density function distribution of the natural streamflow reconstruction performance expressed using CC, NSE, and KGE skill metrics (all range from 0 to 1) for monthly streamflow in 425 catchments throughout mainland China. The CC values are generally high, ranging from 0.46 to 0.99 with a mean of 0.90 for the 425 catchments (Figure 2d), indicating that the model properly captures the trends in natural streamflow dynamics. For the NSE and KGE metrics, respectively, about 77% and 92% of the 425 catchments have NSE and KGE skill scores >0.7. Overall, these skill metric value distributions confirm that our study provides a reliable reconstruction of the natural streamflow for these 425 catchments across China.
The climate variables, including precipitation, temperature, and potential evapotranspiration were calculated and then aggregated based on the climate data sets and the catchment boundaries ( Figure S2 in Supporting Information S1). We also inspected the long-term trends of temperature and precipitation for the 425 catchments across the 10 river basins. Linear trend analysis of annual climate variable anomalies shows a significant warming (median trend = 0.2°C/10 yr) but an insignificant wetting trend (median trend = 0.4%/10 yr) for the 425 catchments during the period of 1961-2018 ( Figures S3 and S4 in Supporting Information S1), which is consistent with previous climate research results, such as Xi et al. The precipitation and temperature changes in warm and cool seasons across the 425 catchments are shown in Figure S4 in Supporting Information S1. Results show that both precipitation and temperature have asymmetrical seasonal changes, with higher warming and wetting trends in cool seasons than in warm seasons ( Figure S4 in Supporting Information S1), which is consistent with a previous study by Piao et al. (2010), who reported a rate of winter warming about four times higher than the rate of summer warming. At basin scale, the trends of several climate variables varied across the 10 river basins (Figure 3). Water year precipitation of the Hai River basin and the Liao River basin during 1961-2018 generally decreased, with median rates of −1.4%/10 yr and −1.0%/10 yr, respectively, while it greatly increased in the Southeast and Northwest River drainage systems, with median values of 3.3%/10 yr and 2.3%/10 yr, respectively (Figure 3a). Temperature increased in all basins, but with different amplitudes of the warming trend (Figure 3b). The annual warming trend was highest in the Yellow River basin (median trend = 0.3%/10 yr) and lowest in the Yangtze River basin (median trend = 0.1%/10 yr). Similarly, the basin-scale seasonal changes in climate variables are consistent with changes at the national scale-that is, there are more significant climate changes in cold seasons than in warm seasons ( Figure 3). Overall, climate dynamics show large differences across different basins and seasons; thus it is necessary to investigate how streamflow responds to seasonal climate change across different hydroclimatic zones.

Streamflow Sensitivity to Climate Change
Based on reconstructed and calculated long-term natural streamflow and climatic variables over the 425 catchments, we applied an empirical linear statistical model to estimate the streamflow sensitivity coefficient to climate change annually and for warm and cool seasons. The values ε P,a (P a , Q a ), ε P,w (P w , Q w ), and ε P,c (P c , Q c ) indicate, respectively, the annual, warm-season, and cool-season streamflow sensitivity to precipitation change; ε T,a , ε T,w , and ε T,c are the same, except they refer to temperature change. Tperf, Pperf, Tsub, and Psub are intended to further quantify the relative strength of streamflow sensitivity to climate change in warm versus cool seasons (Section 2.4). Figure 4 shows the sensitivity of annual streamflow responses to annual precipitation (ε P,a ) and temperature (ε T,a ), and Figure 5 shows the sensitivity results at basin scale. As an example, an ε P,a of 1.0 in Figure 2a implies that a 1-unit increase in normalized annual precipitation throughout the year resulted in a 1-unit increase in normalized annual streamflow, ε T,a .
At the annual scale, we can see that streamflow sensitivity to precipitation change (mean ε P,a = 0.92, P < 0.1 for all 425 catchments) is much greater than streamflow sensitivity to temperature change (absolute mean ε T,a = 0.05) (Figure 4), confirming the dominant effect of precipitation on land surface processes (Berghuijs et al., 2017;Risbey & Entekhabi, 1996;Syed et al., 2004). Consistent with this, Yang and Yang (2011) reported precipitation elasticity 3-10 times greater than the other climatic variable elasticities in the Hai River and the Yellow River basins of China. The ε P,a and ε T,a coefficients have spatial differences with magnitudes ranging from 0.41 to 1.32 and from −0.54 to 0.38, respectively. The values for annual streamflow sensitivity to precipitation change, ε P,a , are all positive for the 425 catchments, with the highest in the Southeast River drainage system (median ε P,a = 1.0, P < 0.1) and the lowest in the Yellow River basin (median ε P,a = 0.85, P < 0.1), indicating that annual streamflow increases or decreases as annual precipitation increases or decreases (Figures 4a and 5a). But streamflow sensitivity to temperature change shows different response directions-some positive and some negative-across the 425 catchments (Figure 4b). For example, annual streamflow for most catchments of the Songhua River basin shows significant positive sensitivity to changes in temperature (median ε T,a = 0.06, Figures 4b and 5b), that is, streamflow increases under a warming climate, while significant negative sensitivity is apparent in the Yellow River basin (median ε T,a = −0.02, Figures 4b and 5b), that is, streamflow decreases with warming. Figure 6 shows the seasonal streamflow sensitivity to precipitation change. The sensitivity shows large differences between warm and cool seasons; namely, sensitivity in warm seasons (mean ε P,w = 0.90) is significantly greater than in cool seasons (mean ε P,c = 0.46), and in cool seasons, negative precipitation sensitivity can even been seen in certain catchments (Figures 6a and 6b). A previous study noted that the seasonality of streamflow sensitivity to precipitation change actually reflects the time of year when precipitation is greatest (Vano et al., 2015). In our study, water year ε P,a is primarily determined by the variation of precipitation in the warm season ( Figure 5a and Figure S5 in Supporting Information S1), because the East Asian summer monsoon prevails in China (Miao et al., 2019;Zhang et al., 2010), and thus precipitation is generally concentrated in warm seasons (May-October). Catchments with negative precipitation sensitivity in cool seasons are mainly located in the upper sections of the Yangtze River basin and the Huai River basin (Figure 6b). The reason for this negative sensitivity may be that the cool-season precipitation in these catchments is too low, so the streamflow variation cannot be well explained by the precipitation change. However, most catchments exhibit consistent positive streamflow sensitivity to precipitation changes in the warm and cold seasons (Pperf > 0, purple catchments in Figure 6c). We also found that most catchments across China have exhibit greater streamflow sensitivity to precipitation changes in warm seasons than in cool seasons (indicated by the positive Psub values in Figure 6d), except for some catchments in southeastern regions. Figure 7 shows the streamflow sensitivity to temperature change in warm and cool seasons. Similar to what we see for precipitation, the streamflow sensitivity to temperature change in warm seasons is more consistent with the results of annual temperature sensitivity than is the sensitivity in cold seasons (Figure 4a, Figures 7a and 7b, and Figure S5b in Supporting Information S1). We found that streamflow sensitivity to temperature change in the cold season is considerably greater than in the warm season, particularly in the Northwest River drainage system ( Figure S5b in Supporting Information S1 and Figure 5b). Meanwhile, the streamflow response to seasonal temperature change exhibits opposite sensitivity in the cold and warm seasons in most basins (Tperf < 0, Figure 7) except for the northern Songhua River basin (ε T,w > 0, ε T,c > 0, Tperf > 0, P < 0.1). We further quantify the relative strength of streamflow sensitivity to temperature change for warm versus cold seasons using the Tsub index. Contrary to what we see for precipitation, the sensitivity of streamflow to temperature was stronger in the cool season than in the warm season in many catchments (Tsub < 0, Figure 7d). A possible reason for the temperature response's seasonal asymmetry is that the historical warming rate has been higher in cold seasons than in warm seasons ( Figure S4b in Supporting Information S1), and a greater warming rate has the potential to melt more snow in cool seasons (Barnett et al., 2005), especially in the northern regions of China.

Impact Factors of Streamflow Sensitivity to Climate Change
To assess the potential impact factors of the sensitivity of streamflow to climate change across the 425 catchments, we first calculated the averaged values of six catchment characteristics for each catchment ( Figure S6 in Supporting Information S1), including aridity index (P/PET), runoff efficiency (Q/P), elevation (DEM), vegetation cover (NDVI), field capacity (W, unit: mm), and saturated hydraulic conductivity (Ks, unit: mm/day). We then fitted locally smoothing curves (LOWESS, span parameter set to 1) and conducted correlation analysis to explore the relationship between the streamflow sensitivity index and the surrounding environment in a given catchment. As shown in Figure 8, we found that the water year streamflow sensitivity to precipitation change (ε P,a ) shows significant negative correlations with the aridity index (P/PET, Pearson's r = −0.44, P < 0.01) and saturated hydraulic conductivity (Ks, Pearson's r = −0.37, P < 0.01), while showing positive correlations with the runoff efficiency (Q/P, Pearson's r = 0.18, P < 0.01), vegetation cover (NDVI, Pearson's r = 0.33, P < 0.01), and field capacity (W, Pearson's r = 0.38, P < 0.01). The relationship between catchment characteristics and streamflow sensitivity to changes in temperature is not as strong as the relationship between catchment characteristics and the streamflow sensitivity to changes in precipitation. Streamflow sensitivity to temperature change (ε T,a ) is strongly related to the vegetation conditions and field capacity when the sensitivity index is significantly positive (Figure 8), while irrelevant in catchments with negative ε T,a .
On a seasonal scale, the relationship between warm-season streamflow sensitivity to climate change and catchment characteristics is similar to the results for the water-year scale, but there are significant differences in cool seasons, especially with P/PET, DEM, and NDVI ( Figures S7 and S8 in Supporting Information S1). We found that the magnitude of streamflow sensitivity to temperature change in cool seasons, ε T,w , has a significant positive correlation with DEM (Pearson's r = 0.31, P < 0.01); that is, the higher the elevation, the greater the ε T,w ( Figure S8 in Supporting Information S1). A possible reason for this positive relationship between elevation and streamflow sensitivity to change in cool-season temperature is that the strong cool-season warming ( Figure S4 in Supporting Information S1) has increased the snowmelt and glacial melt volume in the high-elevation catchments where snow and glaciers are prevalent. We also find that streamflow sensitivity to precipitation change in cool seasons is significantly correlated to the climate conditions in a given catchment; that is, a catchment in a warm and humid region has higher cool-season streamflow sensitivity to precipitation change than catchments in cool and dry (water-limited) regions ( Figure S9 in Supporting Information S1), which broadly agrees with the findings of Vano et al. (2015) for the Pacific Northwest.
To explore the streamflow sensitivity of catchments to shifts in seasonality, we defined two indexes, Tsub and Psub (Equation 9 in Section 2.4), to further quantify the relative strength of streamflow sensitivity to climate change in warm versus cool seasons. We plotted the Psub and Tsub indexes across the 425 catchments in relation to six catchment characteristics for those catchments and LOWESS-smoothed (span parameter = 1) results in Figures 9a and 9b.
We found that the Psub index ( Figure 9a) has significant positive correlations with DEM (Pearson's r = 0.43, P < 0.01) and PET/P (Pearson's r = 0.55, P < 0.01), while it has significant negative correlations with the Q/P (Pearson's r = −0.52, P < 0.01), NDVI (Pearson's r = −0.32, P < 0.01), and W (Pearson's r = −0.56, P < 0.01) indexes, indicating that an increase of the Psub index is associated with higher elevations, drier conditions, less Figure 8. Relationships between water year streamflow sensitivity to climate change (ε P,a and ε T,a ) and catchment characteristics, including (a) aridity index (P/PET), (b) runoff efficiency (Q/P), (c) elevation (DEM), (d) vegetation cover (NDVI), (e) field capacity (W, unit: mm), and (f) saturated hydraulic conductivity (Ks, unit: mm/day). Each plot shows LOWESS-smoothed curves of scatterplots (points are not plotted for clarity) between the streamflow sensitivity to climate changes and physical characteristics for the 425 catchments. The LOWESS smoothing span parameter is 1.0, and shading denotes the 95% confidence interval for the possible locations of the smoothed lines. To simplify, subplots (b-f) use the same legend as subplot (a), and it is not drawn separately.
vegetation, and lower infiltration rates in a given catchment. In other words, in drier and higher-altitude catchments, their streamflow sensitivity to precipitation in different seasons is more likely to present asymmetrically and to exhibit lower streamflow sensitivity to precipitation change in cold seasons. This is understandable in the context of climate seasonality: higher, drier, less vegetated, lower-infiltration catchments are dominated by highland and continental climate regimes and thus generally exhibit strong precipitation seasonality. Therefore, uneven seasonal precipitation enhances the asymmetry of streamflow response to seasonal precipitation change, leading to a stronger inclination toward warm-season precipitation changes in those water-limited regions (i.e., more positive Psub). Figure 9b shows the relationship between the Tsub index (asymmetry of seasonal streamflow sensitivity to temperature change) and six catchment characteristics. Results show that the Tsub values are generally equal to zero (absolute sensitive values are the same in warm and cold seasons) in those catchments with low elevation and high vegetation cover (Figure 9b). With an increase of elevation and a decrease of vegetation coverage, the seasonal asymmetry effect of streamflow sensitivity to temperature change becomes more obvious but is more nonlinear; that is, the season during which streamflow sensitivity to temperature change is found to be larger is roughly the same. We also found that the curves of W, Q/P, and PET/P share similar bimodal patterns. Overall, the relationship between the catchment characteristics and the streamflow response asymmetry to seasonal warming is more complex than for precipitation.

Conclusion
Exploration as to how streamflow sensitivity to climate change varies with hydroclimate regime and season is essential to better understand the interactions between land surface hydrology and climate change. In this study, we first used a hydrological model and a post-processing bias-correction algorithm to reconstruct long-term natural streamflow for 425 catchments across China. The reconstructed natural streamflow was used along with precipitation and temperature series for the basin-averaged water year and warm and cool seasons to construct a multivariate linear regression model to analyze the spatial pattern and seasonality of streamflow sensitivity to climate change. Six average catchment characteristics-including aridity index, runoff efficiency, elevation, vegetation cover, field capacity, and saturated hydraulic conductivity-were derived and used to assess the potential impact factors of the sensitivity of streamflow to climate change across the 425 catchments.
Based on our work, we conclude the following: 1. The reconstructed streamflow matched well with the natural streamflow, with Pbias values within the ±10% acceptable range for most catchments (402 of 425 catchments). Climate change analysis shows significant warming (median trend = 0.2°C/10 yr) and insignificant wetting (median trend = 0.4%/10 yr) trends for the 425 catchments during the period of 1961-2018. Both precipitation and temperature have asymmetrical seasonal change patterns with higher change trends in cool seasons than in warm seasons. 2. At an annual scale, the streamflow sensitivity to precipitation change (mean ε P,a = 0.92, P < 0.1 for all 425 catchments) is much larger than the streamflow sensitivity to temperature change (absolute mean ε T,a = 0.05). And the water year streamflow sensitivity to precipitation change values are all positive for the 425 catchments, with the highest in the Southeast River drainage system and lowest in the Yellow River basin. Yet the streamflow sensitivity to changes in temperature shows different positive and negative responses across the 425 catchments, including significant positive sensitivity for most catchments in the Songhua River basin and significant negative sensitivity in the Yellow River basin. 3. Seasonal sensitivity evaluation results indicated that streamflow sensitivity to precipitation change in warm seasons (mean ε P,w = 0.90) is significantly larger than in cool seasons (mean ε P,c = 0.46), and most catchments have a consistent positive streamflow sensitivity to precipitation change in warm and cold seasons. However, the sensitivity of streamflow to changes in temperature is stronger in the cool season than in the warm season, particularly in the Northwest River drainage system. 4. Catchments with higher elevation, drier climate, less vegetation cover, and lower infiltration rates are more prone to show asymmetric seasonal streamflow sensitivity to precipitation change, with lower sensitivity in cold seasons. The magnitude of streamflow sensitivity to temperature change in the winter season, ε T,w , is significantly affected by DEM values; that is, the higher the elevation, the greater the streamflow sensitivity to temperature change in cool seasons.