Quantifying Impacts of Forest Recovery on Water Yield in Two Large Watersheds in the Cold Region of Northeast China

Abstract

In northern China, large-scale reforestations were implemented to restore the ecosystem functions (e.g., hydrology function). However, few studies have been conducted to quantify the relative contributions of forest recovery to water yield in boreal forest region across the globe. In this study, the impacts of forest recovery on the changes in mean annual water yield were assessed in two large forested watersheds in the boreal forest region of northeast China using three different approaches. As commonly considered, the results confirmed that forest recovery was the dominant driver of the reductions in annual water yield in the two watersheds in the past three decades (1987–2016), explaining 64.3% (15.4 mm) and 87.4% (40.7 mm) of variations in annual water yield for Upper Tahe watershed (UTH) and Xinancha watershed (XNC), respectively. By contrast, climate variability played minor role in annual water yield variation, explaining only 35.7% (8.5 mm) and 12.6% (7.2 mm) for UTH and XNC, respectively. The response differences between the two watersheds may mainly be attributed to differences in forest type, topography and climate regimes. This study provided important insight into sustainable forest and water resources management in the region.

1. Introduction

The large-scale reforestation programs have greatly increased forest coverage in China since the 1980s, restoring ecosystem services and benefits of the communities [1]. Despite the positive effects of reforestation in ecosystems restoration, the massive afforestation was also considered a threat to the water resources or supply in some regions [2]. For example, Qiu [3] reported that reforestation may potentially contribute to the droughts in Southwest China. Although numerous existing studies were dedicated to assessing the impacts of reforestation on water yield globally with varying numbers of watersheds and watershed sizes, as reviewed by Zhang et al. [4] and Li et al. [5], there exists limited studies [6,7] examining the effects of afforestation and reforestation on water yield in boreal regions, and, particularly, a lack of studies in the boreal forest region in China [8]. This raises a critical need to study such effects in the boreal regions to enrich our knowledge of the relationship between forest recovery and water resources.

In addition to forest recovery, climate change is also a critical driver to the hydrology in large forested watersheds [9,10,11]. For instance, Li et al. [5] reported that forest and climate played a co-equal role in hydrological variations. Wei et al. [11] found that forest change can only explain 30% of hydrological variations in the global forested regions in the period of 2000–2011. This further highlighted the importance of considering both forest and climate change in watershed assessment.

To quantify the relative contribution of forest cover change to the changes in water yield, the effects of climate variability have to be removed or quantified firstly. Several recent studies have successfully separated the relative impacts of forest cover change and climate variability on water yield in large watersheds using different methods. For instance, Wei and Zhang [9] and Zhang and Wei [12] applied Modified Double-Mass Curve (MDMC) to quantify the relative contributions of the changes in forest cover to water yield in large forested watersheds in Canada and China, respectively. Zhao et al. [13] used Time Trend Analysis (TRA) and Sensitivity-Based Methods (SBM) to separate the effects of vegetation change and the effects of climate variability on streamflow in different catchments in New Zealand, South Africa and Australia, respectively. The Budyko Framework (BF) is also widely used to quantify the impacts of vegetation changes on streamflow [14,15,16]. Wei et al. [17] compared pros and cons of eight methods for quantifying the relative contribution of forest or land cover change and climatic variability to hydrology in large watersheds and concluded that each method must be supplemented by other methods to achieve a robust conclusion. In this study, we applied three methods including MDMC, TRA and SBM to get the reliable results.

The cold region of Northeast China can be typically represented by the boreal coniferous forest in the Da Hinggan Mountains [18] and the boreal/temperate transition mixed coniferous and broadleaved forest in the Xiao Hinggan Mountains [19], respectively. These forests experienced long-term timber harvesting from the 1960s to 1990s and a forest recovery period through natural forest protection projects in the late 1990s [20]. As a result, the forest coverage has greatly increased in the past three decades, which provides a suitable base for investigating the interactions between forest recovery and water yield in this region.

Two large forested watersheds including Upper Tahe (UTH) watershed and Xinancha (XNC) watershed in the cold region of Northeast China were selected for this study. Long-term (1987–2016) forest cover and hydrometeorological data were used: (1) to quantify the relative contributions of forest recovery to the long-term water yield; (2) to examine the sensitivity of annual water yield to forest recovery; and (3) to explore the implications of our research findings for watershed management in the cold region.

2. Materials and Methods

2.1. Study Watersheds

The UTH watershed and XNC watershed have drainage areas of 2359 and 2582 km², respectively, and are located in Heilongjiang Province in the high latitude cold region of Northeast China (Figure 1a). According to Chinese vegetation database administrated by Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC), the UTH watershed is zoned in boreal coniferous forests, while XNC watershed is zoned in the boreal/temperate transition mixed coniferous and broadleaved forest (Figure 1b). The two watersheds are both dominated by gentle hills with the average slope of 12.5° and 13.4°, respectively, and the elevations range from 96 to 1276 m above sea level (Figure 1c,d). Both watersheds are dominated by dark brown earths and brown coniferous forest soils. The main characteristics of two watersheds are summarized in

Table 1. Watershed characteristics for the UTH watershed and XNC watershed.

Metrics	UTH Watershed	XNC Watershed
Drainage area (km2)	2359	2582
Mainstream length (km)	83	82
Average elevation (m)	783	501
Elevation range (m)	432–1276	69–1226
Average slope (°)	12.5	13.4
Soil type	Dark brown earths and brown coniferous forest soil	Dark brown earths and brown coniferous forest soil
Annual mean precipitation (mm)	534.8	706.0
Annual mean PET (mm)	520.5	711.8
Annual mean air temperature (°C)	−2.1	2.4
Annual mean flow (mm)	297.3	298.7
Average forest cover (%)	75.6	87.7
Forest type	Boreal coniferous forest	Mixed coniferous and broadleaved forest
Hydrometric station	Xinlin	Nancha
Climate stations	Xinlin	Nancha, Xiaobai, Nanlie

. According to the Circum-Arctic Map of Permafrost and Ground Ice Conditions, Version 2 (http://nsidc.org/data/ggd318#) [21], the UTH watershed spans discontinuous and sporadic permafrost zones, while the XNC watershed spans seasonal frozen ground without permafrost. It should be noted that warming climate might further enhance permafrost thaw particularly in boreal regions [22,23]. The permafrost thaw can alter the regional hydrological processes [24,25]. Therefore, to minimize the effects of permafrost thawing on streamflow, the period from 1987 to 2016 with a relative stable mean annual air temperature, as found by Duan et al. [8] and Duan et al. [26], was selected to quantify the relative contribution of forest recovery on water yield in this study.

The two study watersheds are characterized by a typical continental monsoon climate. Based on the climate data from 1987 to 2016, the average annual precipitation of two watersheds are 534.8 and 706.0 mm for the UTH and XNC, respectively, of which approximately 85% occurs as rain from May to September (the wet season) (Figure 2). The average annual air temperature for the UTH and XNC watersheds is −2.1 and 2.4 °C, respectively, with the average highest and lowest occurring in July and January, respectively (Figure 2). The UTH watershed is located in the boreal coniferous forest zone, where the native vegetation consists of forest communities dominated by larch (Larix gmelinii), along with broadleaf species, such as birch (Betula platyphylla). The XNC watershed is in the boreal/temperate transition mixed coniferous and broadleaved forest zone, where the flora is more diverse [18]. In this region, the coniferous forest is dominated by Korean pine (Pinus koraiensis), and the broadleaved species mainly include Tilia amurensis, Fraxinus mandschurica, and Mongolian oak (Quercus mongolica). The forests were mainly regenerated from the natural forest succession with the major secondary forests of birch in UTH watershed and Mongolian oak forest in XNC watershed.

2.2. Data

2.2.1. Forest Cover Data

Annual forest cover data from 1987 to 2016 were obtained from the forest resource inventory database administrated by the Xinlin Forestry and Yichun Forestry Bureaus [8,27], respectively. Following the protocol of “Observation Methodology for Long-term Forest Ecosystem Research” of the National Standards of the People’s Republic of China (GB/T 33027-2016), the forest coverage was calculated as percent ratio of all forest area with canopy coverage greater than 30% over the total area of the study watershed. As shown in Figure 3, the forest cover of two study watersheds showed a significantly positive trend (p < 0.001, Mann-Kendall test) during the entire study period. The Pettitt’s Tests [28] indicated that the forest coverage of UTH and XNC watersheds had the statistically significant change points in 2002 and 2000 (p < 0.05), respectively. The forest coverage of UTH watershed did not show very much change from 1987 to 2002 with the average of 70.2%, but increased from 72.2% to 87.5% from 2003 to 2016 (Figure 3). The forest coverage of XNC watershed increased from 75.8% to 86.6% from 1987 to 2000, and showed a slightly increasing trend from 2001 to 2016 with the average of 88.4%.

2.2.2. Hydrometeorological Data

The annual and seasonal streamflow data from 1987 to 2016 were both produced by the aggregation of daily streamflow data, which were calculated based on the relationship between streamflow and discharge stage height collected at the Xinlin and Nancha hydrometric stations (Figure 1), respectively. The details on annual streamflow, Annual peak flows, the monthly precipitation, monthly mean, maximum and minimum air temperatures, and monthly sunshine hours from 1987 to 2016 are shown in Table 1. The averaged watershed-based precipitation estimates were derived by the Thiessen polygon method for XNC watershed. To understand the relative contribution of seasonal streamflow variations to annual streamflow changes, streamflow data from the entire year was summarized into four seasons: spring (April to May), summer (June to August), autumn (September to October), and winter (November to March).

2.3. Data Analysis Methods

2.3.1. Trend Analysis on Hydrological and Climatic Series

Trend analysis is useful for understanding dynamics of hydrological and climatic variables over a long-term period [29,30]. A nonparametric trend test, namely the Mann-Kendall test [31,32], was adopted to detect whether significant trends exist in the long-term hydrological and climatic data [29]. The magnitude of the trend, β, or the slope (change per unit time) was used to describe the change rate of long-term hydrological and climatic variables, which was estimated using the Sen’s slope method [33,34] as follows:

$$\beta =Median\left[\frac{X_j-X_i}{j-i}\right]\mathrm{for}\mathrm{all}i<j$$

(1)

where 1 < i < j < n and is the median of all possible combinations of pairs for the entire data set.

2.3.2. Time Series Correlation Analysis

Annual time series of hydrological and forest cover data were used to conduct cross-correlation analyses for the entire study period and detect statistical relationships between these data sets. To eliminate the autocorrelations in the data series, all data series were pre-whitened using the best fitting autoregressive integrated moving average models (ARIMA) [35,36]. Then, residual series were used in the cross-correlation analysis. These analyses were conducted using the software STATISTICA 7 (StatSoft^®, Palo Alto, CA, United States).

2.3.3. Separation of the Impacts of Climatic Variability and Forest Recovery on Annual Streamflow

Modified Double Mass Curve Method

Modified Double Mass Curve (MDMC) developed by Wei and Zhang [9] was firstly used to separate the relative contributions of climate variability and forest recovery on annual streamflow and detect the break point in annual streamflow caused by forest recovery. The modified double mass curve assumes that there is a linear relation between streamflow and effective precipitation, and plots accumulated annual streamflow versus accumulated annual effective precipitation. The annual effective precipitation refers to the difference between annual precipitation and annual evapotranspiration [9]. In the period when there are no significant impacts of forest cover change on annual streamflow, a straight line should be expected to describe the relation between annual streamflow and annual effective precipitation. However, with a management or other disturbances in the forest, a curve line with a turning point would be expected. The statistical significance of the break points was confirmed by the interrupted ARIMA model.

Once the statistical significance of the break point was confirmed, the whole study period was subsequently divided into reference (before the break point) and forest recovery (after the break point) periods. The difference between the average observed annual streamflow (Q₂) and the average annual streamflow predicted by the baseline in the forest recovery period can be attributed to the effects of forest recovery on annual streamflow (ΔQ_F). The variation of annual streamflow caused by climate change (ΔQ_C) can be determined as:

$$\Delta Q_C=\Delta Q-\Delta Q_F$$

(2)

where ΔQ is the deviation of average annual streamflow between forest recovery period and reference period:

$$\Delta Q_{}=Q_2-Q_1$$

(3)

where Q₁ and Q₂ are the average annual streamflow in the reference and forest recovery periods, respectively.

The relative contributions of forest recovery and climate variability to the changes in mean annual streamflow can be estimated as:

$$R_F=\frac{\left|\Delta Q_F\right|}{\left|\Delta Q_F\right|+\left|\Delta Q_C\right|}\times 100\%$$

(4)

$$R_C=\frac{\left|\Delta Q_C\right|}{\left|\Delta Q_F\right|+\left|\Delta Q_C\right|}\times 100\%$$

(5)

In this study, potential evapotranspiration is estimated using the temperature-based Hamon method [37], which has been shown to provide reasonable potential evapotranspiration for forested regions [8,38,39]:

$$PET=0.1651\times D\times V_d\times K\times N$$

(6)

where D is the monthly average time from sunrise to sunset (from the climate station) in multiples of 12 h, V_d is the saturated vapor density (g m⁻³) at the monthly mean air temperature (T, °C) as shown in Equation (7), and N is the number of days in each month.

$$V_d=216.7\times \frac{V_s}{T+273.3}$$

(7)

V_S is the saturated vapor pressure (millibars), expressed as:

$$V_S=6.108\times \exp \left[17.26939\times \frac{T}{T+273.3}\right]$$

(8)

where K is a correction coefficient to adjust PET from the Hamon’s method to reflect realistic values for PET. The reported K values range from 1.0 to 1.4 [37,39]. Due to the expected low potential evapotranspiration in the cold study watersheds, the K value was set to be 1.1 as suggested by Duan et al. [8] and Sun et al. [40].

In this study, annual evapotranspiration (ET) was firstly estimated using both the Budyko (Equation 8) and Zhang equations (Equation 9) and, then, compared with the difference between the long-term mean annual precipitation and mean annual streamflow, which is considered as the “actual evapotranspiration” because that it is reasonable to assume that changes in soil water storage are zero and the changes in the recharge to groundwater are small over a long period of time (i.e., 5–10 years) in the catchment water balance framework [41]. We found that the annual ET calculated by Budyko equation was closer to the actual evapotranspiration. Thus, the results from Budyko equation was used in this study.

$$ET={\left\{P\left[1-\exp (-PET/P)\right]\times PET\times \tanh (P/PET)\right\}}^{0.5}$$

(9)

$$ET=P[1+w(PET/P)]/[1+w(PET/P)+P/PET]$$

(10)

where w is plant available water coefficient, which is set to be 2 given the large proportion of forest cover in study watersheds [41]. P and PET represent annual precipitation and annual potential evapotranspiration, respectively.

Sensitivity-Based Method

The variations in mean annual streamflow attributed to climate variability (ΔQ_C) is calculated from the changes of annual precipitation (P) and potential evapotranspiration (PET) between the different phases in the sensitivity-based method as follow [42,43]:

$$\Delta Q_C={\psi }_P\Delta P+{\psi }_{PET}\Delta PET$$

(11)

where ΔP and ΔPET are the changes in precipitation P and PET between the forest recovery period and the reference period. ${\psi }_P$ and ${\psi }_{PET}$ are streamflow sensitivity coefficients to P and PET as expressed below [44]:

$${\psi }_P=\frac{1+2x+3wx}{{(1+x+w{x}^2)}^2}$$

(12)

$${\psi }_{PET}=-\frac{1+2wx}{{(1+x+w{x}^2)}^2}$$

(13)

where x is the mean annual index of dryness and is equal to PET/P and w is plant available water coefficient the same as Equation (10).

Once ΔQ_C was estimated, the ΔQ_F can be determined by the equation as follow:

$$\Delta Q_F=\Delta Q-\Delta Q_C$$

(14)

Time Trend Analysis Method

In the time trend analysis method, the Kendall–Theil robust line method [45] was firstly used to create a linear equation between annual streamflow and annual effective precipitation in the reference period as Equation (15), and the equation was then applied in the forest recovery period to predict annual streamflow. Thus, the difference between mean annual streamflow (Q₂) and mean annual predicted streamflow (Q_2P) in the forest recovery period was considered the effects of forest recovery on annual streamflow (ΔQ_F), as shown in Equation (16). Once the ΔQ_F was estimated, the ΔQ_C can be determined by Equation (2).

$$Q_t=a{P}_e{}_t+b$$

(15)

$$\Delta Q_F=Q_2-Q_{2P}$$

(16)

where Q_t and P_et are annual streamflow and annual effective precipitation at the tth year in the reference period, respectively, and a and b are the constants of linear equation.

3. Results

3.1. Trends of Annual and Seasonal Hydrometeorological Variables

For annual hydrometeorological data series from 1987 to 2016 (

Table 2. Results of Mann-Kendall trend tests on hydrometeorological variables in the UTH watershed and XNC watershed from 1987 to 2016.

Period	Watershed	Q		P		T		PET
		Slope 1 (mm/year)	p	Slope (mm/year)	p	Slope (°C/year)	p	Slope (mm/year)	p
Annual	UTH	−2.0	0.26	−1.2	0.48	0.00	0.93	0.58	0.15
	XNC	−1.1	0.67	0.25	0.92	0.00	0.75	0.59 *	0.01
Spring	UTH	0.35	0.55	1.03	0.11	0.02	0.31	0.30	0.11
	XNC	1.09 *	0.02	2.36 *	0.01	0.02	0.33	0.25	0.17
Summer	UTH	−1.88	0.32	−1.90	0.32	0.02	0.16	0.58	0.08
	XNC	−0.04	1.00	0.31	0.97	0.03	0.07	0.47	0.11
Autumn	UTH	−0.29	0.62	−0.16	0.86	−0.02	0.41	0.00	0.94
	XNC	−0.63	0.50	−1.67	0.18	0.01	0.72	0.10	0.27
Winter	UTH	−0.04	0.44	−0.2	0.36	−0.03	0.25	−0.04	0.83
	XNC	−0.03	0.86	1.18 *	0.02	−0.04	0.38	−0.11	0.55

* values indicate the statistical significance at the level of 0.05. ¹ The Slope was estimated using the nonparametric median-based slope method [33,34]. Q, streamflow; P, precipitation; T, mean annual air temperature; PET, potential evapotranspiration ET, evapotranspiration; Spring, April to May; Summer, June to August; Autumn, September to October; Winter, November to March.

), very few significant (p < 0.05) trends were found, mainly on annual PET, spring flow rate, and spring and winter precipitations in XNC watershed in comparison with none significant (p > 0.05) trend in UTH watershed. Although there was no significant trend in annual streamflow in UTH watershed, the annual flow rate decreased by 2.0 mm/year vs. 1.1 mm/year in XNC watersheds. However, the summer streamflow in UTH watershed was characterized by an apparent decreasing pattern (1.88 mm/year) comparing to the decreasing precipitation (1.90 mm/year).

3.2. Cross-Correlations between Forest Cover and Hydrological Variables

Cross-correlation analysis suggests that annual streamflow in UTH and XNC watersheds was significantly negatively correlated with annual forest cover from 1987 to 2016 (

Table 3. Cross-correlations between forest cover and hydrological variables.

Variables	Forest Cover of UTH (1,1,1)				Forest Cover of XNC (1,1,1)
	ARIMA 1 Model	Cross-Correlation Coefficient	p	Lag	ARIMA Model	Cross-Correlation Coefficient	p	Lag
Annual	(1,0,0)	−0.51 **	0.005	9	(1,0,0)	−0.54 **	0.001	5
Spring	(1,0,0)	0.46 *	0.019	4	(0,0,1)	0.41 *	0.026	0
Summer	(0,0,1)	−0.40 *	0.038	10	(1,0,0)	−0.46 **	0.007	5
Autumn	(1,0,1)	−0.41 *	0.017	9	(0,0,1)	−0.40 *	0.025	5
Winter	(1,0,0)	−0.28	0.170	5	(1,0,0)	0.16	0.404	0

¹ ARIMA is autoregressive integrated moving average. ** values indicate the statistical significance at the level of 0.01. * the statistical significance at the level of 0.05.

) except for winter season regardless of watersheds. However, there is a significantly positive correction between the streamflow and forest cover during spring season for both watersheds. In addition, the lags between hydrological variables and forest cover varied with watershed, 4-10 years for UTH and 0–5 years for XNC watersheds, respectively.

3.3. Separating the Relative Contributions of Forest Recovery and Climate Variability to the Changes in Annual Streamflow

Break points were detected in the modified double mass curves (MDMCs) in 2003 and 2001 for UTH and XNC watersheds (Figure 4), respectively. The fitted interrupted ARIMA model of the MDMCs slopes further confirmed that the break points were statistically significant (p < 0.05,

Table 4. Interrupted ARIMA models for slopes of MDMC in Figure 4.

Model Input	Model Structure	Parameter Estimation
		q(1) 1	Ω 2
Slope of MDMC 3 of UTH in Figure 4a	Interrupted ARIMA: (0,0,1), intervention at year 2003	−0.73 (p < 0.05)	0.52 (p < 0.05)
Slope of MDMC of XNC in Figure 4b	Interrupted ARIMA: (0,0,1), intervention at year 2001	−0.76 (p < 0.05)	0.73 (p < 0.05)

¹ q is moving average parameter, ² Ω is intervention parameter for abrupt permanent intervention type, ³ is modified double mass curve.

). Thus, the entire study period (1987–2016) was then identified as two periods: the reference period (1987 to 2002 for UTH and 1987 to 2000 for XNC) and the forest recovery period (2003 to 2016 for UTH and 2001 to 2016 for XNC).

Forest recovery decreased the mean annual streamflow by 16.9 and 43.4 mm from the reference period in UTH and XNC watersheds, respectively, the corresponding relative contributions to the changes in mean annual streamflow being 70.6% and 93.3%, respectively. By contrast, climate variability induced reductions in mean annual streamflow were 7.0 and 3.1 mm, respectively, with corresponding relative contributions being 29.4% and 6.7%, respectively.

The results of the sensitivity-based method (

Table 5. The relative contributions of climate variability and forest recovery calculated by the sensitivity-based method.

Watershed	Sub-Periods	P (mm)	PET (mm)	ΔP (mm)	ΔPET (mm)	β	γ	ΔQ (mm)	ΔQC (mm)	ΔQF (mm)
	1987–2002	540.7	517.3
UTH				−11.3	7.0	0.58	−0.33	−23.9	−8.8	−15.1
	2003–2016	529.5	524.3
	1987–2000	715.9	652.8
XNC				−7.7	7.9	0.60	−0.36	−46.5	−7.5	−39.0
	2001–2016	708.2	660.7

) indicate that climate variability decreased mean annual streamflow by 8.8 and 7.5 mm from the reference period to the forest recovery period in UTH and XNC watersheds, respectively, in comparison with the reduction of precipitation (−11.3 mm for UTH and −7.7 mm for XNC) and the increase of PET (7.0 mm for UTH and 7.9 mm for XNC). Thus, the contributions of forest recovery to the changes in mean annual streamflow were determined as −15.1 and −39.0 mm in UTH and XNC watersheds, respectively (Table 5).

The results of time trend analysis method shown in Figure 5 indicated that the mean annual observed streamflow was 14.1 and 39.8 mm lower than the predicted values in forest recovery period in UTH and XNC watersheds, respectively. Thus, the relative contributions of climate variability were calculated as the difference between the total changes in mean annual streamflow and the changes attributed to forest recovery (Equation 1), and these values were −9.8 and −6.9 mm for UTH and XNC watersheds, respectively.

The total changes in mean annual streamflow (ΔQ) and the associated components (ΔQ_C and ΔQ_F) calculated by three independent methods in two watersheds are summarized in

Table 6. Long-term annual streamflow changes and their components in UTH and XNC watersheds.

Method	UTH (from 1987–2002 to 2003–2016)					XNC (from 1987–2000 to 2001–2016)
	Changes in Streamflow (mm)			Relative contributions (%)		Changes in Streamflow (mm)			Relative Contributions (%)
	ΔQ	ΔQC	ΔQF	Climate	Forest	ΔQ	ΔQC	ΔQF	Climate	Forest
MDMC 1	−23.9	−7.0	−16.9	29.4	70.6	−46.5	−3.1	−43.4	6.7	93.3
TRA 2		−9.8	−14.1	40.9	59.1		−6.9	−39.6	14.9	85.1
SBM 3		−8.8	−15.1	36.8	63.2		−7.5	−39.0	16.1	83.9
Average		−8.5	−15.4	35.7	64.3		−7.2	−40.7	12.6	87.4

¹ MDMC is Modified Double Mass Curve method, ² TRA is Time Trend Analysis method and ³ SBM is Sensitivity-based Method.

. The results from three methods produced similar results. The relative contributions of forest recovery to the reductions of mean annual streamflow were −15.4 mm (64.3%) and −40.7 mm (87.4%) in UTH and XNC watersheds, respectively, while the relative contributions of climate variability were −8.5 mm (35.7%) and −7.2 mm (12.6%), respectively. Overall, the impacts of forest recovery on long-term annual streamflow variations were much higher than those from climate variability, which indicated that the variations of water yield in the UTH and XNC watersheds were mainly controlled by forest recovery in the past three decades.

4. Discussion

4.1. The Effects of Forest Recovery on Water Yield

Our results indicated that the forest recovery was the dominant driver of the reduction in mean annual water yield, −15.4 and −40.7 mm for UTH and XNC watersheds, respectively. The negative effects of reforestation on water yield are consistent with many study findings in other regions. For instance, Liu et al. [46] found that the reforestation in the Meijiang watershed (6983 km²) covered by subtropical evergreen broad-leaved forest caused an average annual streamflow reduction of 51 mm in the reforestation period (1985–2006). Tuteja et al. [47] found that annual runoff reductions from pine plantations ranged from 22 to 52 mm/year in the different subcatchments of a large catchment in southeastern Australia. The negative effects of reforestation on water yield are possibly due to the increase in evapotranspiration resulting from vegetation development (e.g., increasing leaf areas and root systems) [46]. Forest recovery can also improve soil conditions and enhance the water storage in aquifers [48]. Such changes in hydrological processes consequently result in the reduction of water yield observed from stream. Such reductions in annual water yield caused by reforestation were also reviewed on 73 watershed studies across the globe by Li et al. [5].

Although the forest recovery in this study showed similar impacts on water yield to the previous studies, the sensitivity of annual water yield to forest cover change in these study watersheds is much greater than that in other studies. Our results indicated that 1% forest cover increase resulted in 0.7% of reduction in annual water yield in the boreal coniferous forest watershed (UTH), while 1.8% of reduction in the mixed coniferous and broadleaved forest watershed (XNC). After examining 61 large (>1000 km²) watersheds across the world, Zhang et al. [4] found that in large mixed forest dominated watersheds, 1% forest cover change can result in 0.80% change in annual runoff vs. 0.24% for large coniferous forest dominated watersheds. The difference in the sensitivity of annual water yield response to forest cover change between our study and the study of Zhang et al. [4] may be due to the differences in climatic regimes. The watersheds in this study are located in the high latitude cold region of Northeast China, where the climate is characterized by a typical continental monsoon climate. Approximately 85% of precipitation occurring as rain during the growing season may have been evapotranspired and consequently enhance the negative effects of forest recovery on water yield, as reported by Zhang, Dawes and Walker [41]. This conclusion can be also reinforced by the fact that the total decreases in streamflow during growing season (from June to October) accounted for 84.4% and 77.9% of the total reductions of annual water yield in UTH and XNC watersheds, respectively. Thus, it is safe to say that the impacts of forest recovery on long-term water yield are closely related to its impacts during wet seasons (i.e., summer and autumn in study watersheds), a good consistence with other summer dominant rainfall regions as well [49,50,51]. In addition, the different contribution rates of forest recovery to streamflow change between two watersheds in our study may also be partially related to the water condition difference (Table 1).

4.2. The Effects of Forest Type and Topography on the Response Intensity of Water Yield to Forest Recovery

Our results indicated that 1% forest cover increase can result in 1.8% reduction in annual streamflow in the mixed coniferous and broadleaved forest watershed (XNC) and 0.7% in the boreal coniferous forest watershed (UTH). The distinct response between two watersheds may be because of the difference in forest types, topography and climate regimes. Firstly, forest type can significantly affect the hydrological response to forest recovery [52]. For instance, in the Puget Sound basin covered by mixed coniferous and deciduous forest in northwestern Washington, USA, the annual streamflow response to forest recovery was estimated as 6.8% [53]. However, in the large boreal coniferous forest watersheds in Northeastern Ontario, Canada, there was no definitive changes in annual water yield with forest cover ranging from 8.6% to 25.2% [54]. Meanwhile, Zhang et al. [4] found that the response intensity of annual streamflow to forest cover change in large mixed forest watersheds was approximately three time of that in large coniferous forest watersheds based on the data from 61 large watersheds across the world. These previous studies and ours all suggested that large boreal coniferous forest watersheds may have the relative higher hydrological elasticity in response to forest cover change than the large mixed forest watersheds. In this study, relatively lower transpiration of boreal coniferous forest than the broadleaved forest in the cold regions [55,56] may lead to smaller changes of evapotranspiration in response to forest recovery. In addition, the different mean annual rainfall between two sites in this study can have a strong influence on the change in annual water yield with forest cover change [41,57]. Bosch and Hewlett [58] and Farley, Jobbágy and Jackson [52] found that vegetation change has the largest absolute impacts on water yield in high-rainfall areas. This is confirmed by our findings that the XNC watershed has a higher magnitude of rainfall than UTH watershed, therefore, a greater response intensity of annual streamflow to forest recovery as well. Additionally, the boreal forest watershed (UTH) has a shorter and cooler growing season than the mixed forest watershed (XNC) because of the higher latitude, which consequently results in a lower evapotranspiration in the watershed. All these different characteristics of forest types and interactions with climate in boreal forest watershed attenuate the strength of the impacts of forest cover change on long-term water yield. Nevertheless, more case and modeling studies can definitely help to examine how forest types affect water yield in large forested watersheds.

Topography also plays a critical role in determining hydrological responses to reforestation [59]. As shown in

Table 7. Averaged slopes in two studied watersheds.

Watershed	Percentage of Watershed Area (%)
	Slope >40°	20–40°	15–20°	10–15°	5–10°	<5°
UTH	0.7	22.0	19.6	30.4	3.8	23.4
XNC	0.5	33.4	34.0	4.6	4.6	23.0

, the UTH watershed is characterized by a more gently topography than XNC watershed. XNC watershed has 68% of the watershed area covered by the slopes greater than 15° vs. 44% in UTH watershed, which may have made the UTH a greater hydrological elasticity in response to forest recovery. Although few studies have been dedicated to detect the impacts of different topography on annual water yield, the studies on the effects of topography on specific flow variable [60] can help understand how different topography affects annual water yield. Liu et al. [61] found that hydrological recovery is limited and slower with reforestation in the steeper watershed. Li et al. [62] also found that the topography indices including perimeter, slope length factor, surface area, openness, and topographic characteristic index can be responsive to the streamflow change, in particular, to the low flow variables in snow-dominated regions in the Southern Interior of British Columbia, Canada. This may be because a watershed with gentler topography would likely have a higher water retention ability due to longer flow paths and residence time and consequently enhance the hydrological elasticity [63,64]. Thus, the different response intensity of annual water yield to forest recovery between two study watersheds was partly explained by the hydrological elasticity differentiated by topography in the two watersheds.

Our results of cross-correlations between forest cover and annual flow indicated that there were nine- and five-year lags between forest cover and annual water yield in the UTH and XNC watersheds, respectively (Table 3). Such lagging effects were mainly due to the delayed hydrological responses to forest recovery, because forest recovery may take years or decades to reach a new hydrological equilibrium [65], particularly in the boreal forest that takes much longer time to recover after disturbance than other forest ecosystem in warmer regions [66]. On the other hand, the longer lags in the boreal coniferous forest watershed (UTH) could also partly be due to the stronger hydrological elasticity discussed above.

4.3. The Relative Contributions of Forest Recovery and Climate Variability to Water Yield Variations

Our results from three independent methods indicated that forest recovery was the dominate driver with the relative contributions to the changes in water yield being 64.3% and 87.4% in the UTH and XNC watersheds, respectively. Consistent results were also found in other regions. For instance, Liu et al. [67] studied in the middle and lower reaches of the Yellow River, China, vegetation changes were the main cause for triggering annual runoff changes, which account for approximately 80% from baseline period to changeable period. Similar findings were also reported for the headwaters of the Yellow River basin by Zheng et al. [68]. By contrast, many more studies conclude climate variability has a similar or greater strength of impacts on water yield compared to forest cover change. For instance, the equal hydrological impacts of forest cover change and climate variability were found in large forested watersheds in the central interior of British Columbia, Canada [9], in the the Upper Minjiang River of Yangtze River basin, China [69] and in the upper reach of the Poyang Lake basin, China [46]. Li et al. [10] found that the relative contribution of forest disturbance was only 27% in the Upper Similkameen River watershed situated between Canada and the USA. Similarly, Shi et al. [70] found that the streamflow was more sensitive to climate variability than land cover change in the Upstream of Huai River, China. The relative contributions of forest cover change and climate variability are largely dependent on their magnitude and the characteristics of watersheds. More case studies would help explore how forest cover change and climate variability interactively affect hydrology in large watersheds.

4.4. Implications and Uncertainty

Boreal coniferous forest and boreal/temperate transition mixed coniferous and broadleaved forest are both located in the remote cold region of Northeast China, where the natural forest experienced prolonged timber harvesting during the 1960s to late 1990s before Natural Forest Protection Project was carried out, and experienced a recovery period in the past three decades. Since 2016, timber harvesting has been completely banned in the natural forest. Although such forest management policies have made forest ecosystems well restored, the potential changes in water resources shrink caused by the increasing forest cover had been ignored. Our results indicated that the mean annual water yield decreased by 8.0% and 15.6% due to the increases of 11.6% and 8.7% in forest cover in UTH and XNC watersheds, respectively, in the past three decades. Such great sensitivity of annual water yield to the increases in forest cover should be a great concern, especially in the region covered by mixed coniferous and broadleaved forest, as the downstream of this region is one of the most important crop-producing areas in China. The reduction of annual water yield in the streams may cause a shortage of water for irrigation in downstream. Thus, the policies of forest management should meet the water level that can maintain aquatic functions and ecosystem integrity as well as ensure the water supply for irrigation in downstream area, which need more future quantitative studies to provide more information. For the boreal forest watershed, the reduction of water yield may not be a serious issue due to the relatively low demands of water supply. Nevertheless, the reduction in water flow caused by forest cover change require further investigation, because they are critical for maintaining the dynamics of in-channel and floodplain habitats that play a critical role in sustaining native biodiversity and ecosystem integrity in rivers [71,72], which can be more important than water supply in the boreal forest region in China.

There are several uncertainties in this study. Firstly, although three independent methods were used to quantify the relative contribution of forest recovery to annual water yield and achieved relatively consistent results, the variations of hydrological processes over study period were rarely understood. For example, annual effective precipitation was used to minimize the impacts of annual precipitation on streamflow in the methods of MDMC and time trend analysis method. However, the changes in intra-annual [73] or seasonal climate patterns [74] can significantly affect inter-annual water yield. In addition, reforestation can potentially more or less increase regional precipitation and water availability, and consequently compensate water loss by increased forest evapotranspiration [1,75]. However, these meteorological and hydrological variations were not considered in our quantitative methods. Second, the PET values were estimated by the Hamon method, which is a widely used temperature-based method. However, relative humidity, wind speed, and solar radiation could be affected by reforestation and, in turn, affect the PET [76,77,78]. These factors were not included in the PET estimations in this study. Thirdly, permafrost thaw and seasonal frost changes caused by climate warming can happen even in a period with the relative stable temperature [79,80], and consequently affect the long-term regional water yield [8,25,81]. Although the period 1987 to 2016 had a relative stable temperature, many studies demonstrated that there was a significant warming trend in Northeast China in the past half century [82,83], thus climate change consideration would be necessary for a longer study. In particular, when permafrost warming or frozen ground degradation has already been observed in this region [84,85], which may completely change the response pattern of streamflow to forest recover. However, the impacts of permafrost thaw and frozen ground degradation on water yield were not considered in this study, which need future more process-based studies to investigate.

5. Conclusions

Based on data from two monitored large watersheds, this study proved that forest recovery was the dominant driver to the reduction of mean annual water yield, while the impacts of climate variability were relatively low in the two large forested watersheds in cold region of Northeast China during the past three decades. The relative contributions of forest recovery to the reductions in mean annual water yield were 64.3% (15.4 mm) and 87.4% (40.7 mm) in UTH and XNC watersheds, respectively, while the rest of the reductions in mean annual water yield were attributed to climate variability. We also found that the response intensity of annual water yield response to increasing forest cover in mixed coniferous and broadleaved forest watershed (XNC) was much greater than that in boreal coniferous forest watershed (UTH). It is well known that forest can conserve water and soil resources, therefore, reduce streamflow. However, the reduction of streamflow responding to the increasing vegetation recover may pose an additional issue to the downstream water supply for irrigating agricultural land. A proper trade-off between forest resource protection and proper downstream irrigation water supply must be sought in the future for an effective ecosystem management. These findings are of great importance for both water resource and forest management in large forested watersheds in Northeast China and similar watersheds in other cold regions.