**1. Introduction**

Global climate change and human activity significantly affect the hydrological cycle in cold regions [1–4]. Trends in runoff and the factors that influence them regionally include runoff trends that are positive in certain cold areas. Runoff may be influenced by glacier melt in response to increasing temperature. This phenomenon occurs in the Upper Khovd River Basin of Central Asia [5,6]. Runoff may also be influenced by increasing precipitation. This process occurs in most of Russia and on the south slope of the Altai Mountains in Northwestern China [5,7]. Runoff may also be influenced by the increases in precipitation, glacial meltwater, and permanently frozen soil meltwater that occur with rising temperatures. These effects are observed in the Yangtze River source region, the Nagqu River Basin in the southern part of the Qinghai-Tibet Plateau, the north and south

**Citation:** Liu, S.; Zhou, Z.; Liu, J.; Li, J.; Wang, P.; Li, C.; Xie, X.; Jia, Y.; Wang, H. Analysis of the Runoff Component Variation Mechanisms in the Cold Region of Northeastern China under Climate Change. *Water* **2022**, *14*, 3170. https://doi.org/ 10.3390/w14193170

Academic Editor: Alexander Shiklomanov

Received: 11 September 2022 Accepted: 5 October 2022 Published: 8 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

slopes of the Tianshan Mountains, and the north slope of Qilian Mountain in Northwestern China [4,5,8–10].

Runoff has been decreasing in response to climate change and human activity in the alpine mountains of Northern Eurasia, Central Asia, South Africa, South America, and elsewhere. Runoff has also decreased as a result of high water consumption in areas such as the inland rivers downstream of the Buqtyrma River Basin in Central Asia [11] and the midstream region of the Heihe River Basin in the arid inland river basins of Northwestern China [12]. Runoff into the rivers of cold regions has been decreasing in response to climate change. Some of these rivers were affected mainly by precipitation. They include the rivers of Mount Kilimanjaro in Tanzania, the tropical Andes of South America, and the south slope of the Altai Mountains. In these cases, the proportions of runoff recharge from glacial meltwater are small [5,13–17]. Certain rivers are affected mainly by temperature, such as those in the Northern Rocky Mountains, where flow attenuation may be caused by a reduction in snowpack accumulation at lower altitudes [18,19]. Certain rivers are affected by both temperature and precipitation, such as those in Northwestern China, the northern part of the Qinghai-Tibet Plateau, and the Songhua River Basin (SRB) in Northeastern China [4,5,20,21]. In regions with few glaciers and minimal rainfall, glacial meltwater decreases and evaporation increases with rising temperature. This phenomenon occurs in the headwater region of the Manas River, the north slopes of the Qilian and Kunlun Mountains, the south slope of the Tianshan Mountains in Northwestern China, and the source regions of the Yellow River Basin in the northern part of the Qinghai-Tibet Plateau [4,5,22,23].

Earlier studies [11,12] have investigated variations in runoff in response to climate change and water use in cold regions. Nevertheless, the mechanisms of runoff component variation are poorly understood. The present study used the Songhua River Basin (SRB) in Northeastern China as an example. We analyzed the mechanisms of runoff component variation during the annual freezing, thawing, and non-freeze-thaw periods of the year based on simulations of the soil freeze-thaw and water cycle processes.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The SRB is located in Northeastern China between 41◦420–51◦380 N and 119◦520–132◦310 E (Figure 1). Its elevation is in the range of ~50–2700 m. The SRB covers 557,000 km<sup>2</sup> and spans Heilongjiang, Jilin, Liaoning, and Inner Mongolia Provinces. The Songhua River is the main tributary of Heilongjiang, and it has two sources to the north and the south. The Nenjiang River to the north originates from Yilehuli Mountain, which is a branch of the Daxing'an Mountains. The second source to the south originates from Tianchi Lake in Changbai Mountain. Precipitation dominates the hydroclimatologic regime of the area, and there is no glaciation. The SRB is characterized by seasonally frozen soil. The maximum freezing depth is >200 cm in the basin. The longest freezing period is from October to July of the following year [24]. The average annual precipitation, temperature, and runoff in the SRB were 533.18 mm, 2.95 ◦C, and 629.9 billion m<sup>3</sup> , respectively, between 1956 and 2018.

### *2.2. Data Collection*

Several parameters could be directly measured for each basin. (a) The digital elevation method (DEM) was implemented at an accuracy of 30 m. (b) Daily precipitation, temperature, humidity, wind speed, and sunshine hours were compiled for 51 national meteorological stations in the Songhua River Basin. These meteorological data were released by the National Meteorological Information Center (2019). (c) Land use data for 1990, 2000, and 2005 were measured at 30-m resolution and provided by the Institute of Geography of the Chinese Academy of Sciences, Xinjiang, China. (d) Soils and their characteristic properties were derived from the Second National Soil Census (NSCO, 1979). (e) Water use data for 2000–2018 were acquired from the SRB Water Resources Bulletin (http://www.slwr.gov.cn/, (10 July 2022)). Water use data for 1956–1999 were extrapolated

from the data for population (1956–2018), irrigated area (1956–2018), gross domestic product (GDP) (1956–2018), and water use (1980–2018). (f) Population, irrigated area, and GDP data for the SRB (1956–2018) were obtained from the Statistical Yearbooks of Heilongjiang, Jilin, Liaoning, and Inner Mongolia Provinces (http://tjj.hlj.gov.cn/tjsj, (10 July 2022), http://tjj.jl.gov.cn/tjsj/tjnj/, (10 July 2022), http://tjj.ln.gov.cn/tjsj/sjcx/ndsj, (10 July 2022), http://tj.nmg.gov.cn/tjyw/jpsj/, (10 July 2022)). (g) Groundwater resources data for 1980–2018 were acquired from the SRB Water Resources Bulletin (http://www.slwr.gov.cn/, (10 July 2022)). *Water* **2022**, *14*, 3170 3 of 15

**Figure 1.** Overview of the SRB. **Figure 1.** Overview of the SRB.

*2.2. Data Collection* The model was verified using measured flow data for the Jiangqiao, Fuyu (1956–2000, 2006–2018), and Jiamusi (1956–2018) Stations.

Several parameters could be directly measured for each basin. (a) The digital eleva-

#### tion method (DEM) was implemented at an accuracy of 30 m. (b) Daily precipitation, tem-*2.3. Research Method*

*2.3. Research Method*

perature, humidity, wind speed, and sunshine hours were compiled for 51 national meteorological stations in the Songhua River Basin. These meteorological data were released by the National Meteorological Information Center (2019). (c) Land use data for 1990, 2000, and 2005 were measured at 30-m resolution and provided by the Institute of Geography of the Chinese Academy of Sciences, Xinjiang, China. (d) Soils and their characteristic The Mann–Kendall trend and Pettitt mutation analyses [25,26] were used to identifying the trends and detect abrupt changes in the measured runoff. The distributed dualistic water cycle model WEP-N was used to simulate the water cycle process in the SRB. WEP-N is driven by multiple factors. A multifactor attribution analysis [27] was used to determine water use and climate change contribution rates to runoff variation in the SRB.

properties were derived from the Second National Soil Census (NSCO, 1979). (e) Water

#### use data for 2000–2018 were acquired from the SRB Water Resources Bulletin 2.3.1. Principles of the WEP-N Model Hydrological Cycle

(http://www.slwr.gov.cn/, (10 July 2022)). Water use data for 1956–1999 were extrapolated from the data for population (1956–2018), irrigated area (1956–2018), gross domestic product (GDP) (1956–2018), and water use (1980–2018). (f) Population, irrigated area, and GDP data for the SRB (1956–2018) were obtained from the Statistical Yearbooks of Heilongjiang, Jilin, Liaoning, and Inner Mongolia Provinces (http://tjj.hlj.gov.cn/tjsj, (10 July 2022), http://tjj.jl.gov.cn/tjsj/tjnj/, (10 July 2022), http://tjj.ln.gov.cn/tjsj/sjcx/ndsj, (10 July 2022), http://tj.nmg.gov.cn/tjyw/jpsj/,(10 July 2022)). (g) Groundwater resources data for 1980– The Water and Energy Transfer Processes and Nitrogen Cycle Processes Model in cold regions (WEP-N) [28] was developed based on the distributed hydrological model in cold regions known as WEP-COR, which couples simulations of natural hydrological and water use processes [24]. The model considered the influences of meteorology, underlying surfaces, and human activity on the water cycle process. The social water cycle simulation included water storage, intake, delivery, use, consumption, and drainage. Water intake and drainage connect the natural and social water cycles.

2018 were acquired from the SRB Water Resources Bulletin (http://www.slwr.gov.cn/,(10 July 2022)). The model was verified using measured flow data for the Jiangqiao, Fuyu (1956– 2000, 2006–2018), and Jiamusi (1956–2018) Stations. The model was calculated using the contour bands inside sub-watersheds. Each contour band was divided according to the land use of the underlying surface, namely, water body, impervious area, soil vegetation, irrigated farmland, or non-irrigated farmland [24,29]. The runoff was calculated based on the proportions of the underlying surface area for each

fying the trends and detect abrupt changes in the measured runoff. The distributed dualistic water cycle model WEP-N was used to simulate the water cycle process in the SRB. WEP-N is driven by multiple factors. A multifactor attribution analysis [27] was used to determine water use and climate change contribution rates to runoff variation in the SRB.

group. Each unit included seven vertical layers. From top to bottom, they were: vegetation canopy or building interception layer, surface depression storage layer, root zone comprised of three layers, transition zone layer, and groundwater layer. The soil below the ground was divided into 11 layers for the soil water and heat coupling calculations. Simulation of the natural water cycle process in the watershed included runoff generation and confluence and comprised the water cycle processes in the surface soil, soil, groundwater, slope ditch, and river channel. The water cycle in the surface soil included precipitation, snow melting, evaporation, infiltration, and surface runoff. The snow melting process was calculated by the temperature-index method [30]. The evapotranspiration values of the surface soil, soil, and water were calculated using the Penman formula [31]. Evaporation from the vegetation canopy was calculated using the Penman–Monteith formula [31]. Infiltration was calculated on the basis of the rainfall intensity and using the Green–Ampt model or the Richards equation [32]. The surface runoff was calculated by the Hortonian or saturation overland flow theory when the precipitation intensity exceeded the infiltration capacity [33]. The soil hydrothermal cycle included heat and water transport during the freezing and thawing periods. Soil heat transport was determined using the basic one-dimensional vertical heat flow movement equation (Equation (1)). Soil water transport was determined using the one-dimensional vertical water flow equation (Equation (2)). The movement of liquid water in the soil was driven by the soil water potential, including pressure, gravity, temperature, and solute potentials (Equations (4)–(6)). The relationship between water and heat transport in frozen soil was characterized as the dynamic balance between the moisture content of the unfrozen water and the negative temperature of the soil. Groundwater movement was calculated with the Boussinesq equation [34]. The exchange between river and groundwater was calculated by Darcy's law and was based on the differences in the water level and the characteristics of the riverbed material [35]. Soil freezing would hinder the exchange between the groundwater and river channels when the temperature exceeds a certain critical value during the freeze-thaw period. The overland confluence was calculated by the kinematic wave method and from the uppermost to the lowermost contour zone of each sub-watershed [36]. Each river channel confluence was calculated using one-dimensional motion waves from upstream to downstream according to the elevation, slope, and Manning roughness coefficient of each contour zone [24]. The influences of temperature on ice formation and river melting were considered.

The temperature difference between the atmosphere and the surface was the heat conduction source. The surface temperature was determined, and the heat flux and temperature of each layer were calculated as follows [37,38]:

$$\frac{\partial}{\partial z} \left[ \lambda\_s \frac{\partial T}{\partial z} \right] = \mathcal{C}\_v \frac{\partial T}{\partial t} - L\_i \rho\_i \frac{\partial \theta\_i}{\partial t} \tag{1}$$

where *z* is the soil layer thickness (m), *T* is the temperature of each soil layer (◦C), *λ<sup>s</sup>* and *C<sup>v</sup>* are the soil thermal conductivity (W/[m· ◦C]) and volumetric heat capacity (J/[m<sup>3</sup> · ◦C]), respectively, *t* is the time (s), and *L<sup>i</sup>* , *ρ<sup>i</sup>* , and *θ<sup>i</sup>* are the latent heat of ice melting (3.35 <sup>×</sup> <sup>10</sup><sup>5</sup> J/kg), the ice density (920 kg/m<sup>3</sup> ), and the volumetric ice content (m3/m<sup>3</sup> ) in the soil, respectively.

The soil liquid water movement was calculated using the Richards equation [39].

$$\frac{\partial \theta}{\partial t} = -\frac{\rho\_l}{\rho\_l} \frac{\partial \theta\_l}{\partial t} - \frac{\partial}{\partial z} \left[ -K(\theta) \frac{\partial H}{\partial z} \right] \tag{2}$$

where *ρ<sup>l</sup>* and *θ* are the soil density (kg/m<sup>3</sup> ) and volumetric water content (m3/m<sup>3</sup> ), respectively, *K*(*θ*) is hydraulic conductivity of the unsaturated soil (m/s), and *H* is the soil water potential (m).

Temperature drives the changes in the water phase. Hydrothermal coupling was calculated as follows:

$$
\theta\_l = \theta\_m(t) \tag{3}
$$

where *θm*(*t*) is the maximum liquid water content corresponding to a negative soil temperature. The soil water potential *H* was calculated as follows [40]:

$$H = z + h\_m + h\_s + h\_t \tag{4}$$

where *h<sup>m</sup>* is the pressure potential (m), *h<sup>s</sup>* is the solute potential (m), and *h<sup>t</sup>* is the temperature potential (m).

The solute potential is a function of temperature and was calculated as follows [40]:

$$h\_s = -\frac{cR}{\mu}T\_K\tag{5}$$

where *c* is the mass of a solute in the unit volume of the solution (kg/m<sup>3</sup> ), *R* is the molar gas constant (8.3145 J/[mol·k]), *µ* is the molar mass of the solute (g/mol), and *T<sup>k</sup>* is the thermodynamic temperature (*K*).

The temperature potential was calculated as follows [41]:

$$h\_t = h\_m G\_{WT} \frac{1}{\gamma\_0} \frac{d\gamma}{dt} \tag{6}$$

where *GWT* is the gain factor (dimensionless), *γ<sup>0</sup>* (*γ<sup>0</sup>* = 71.89 g/s<sup>2</sup> ) is the surface tension of the soil at 25 ◦C (g/s<sup>2</sup> ), and γ is the surface tension of the soil (g/s<sup>2</sup> ). Note that *γ* = 75.6 − 0.1425*T<sup>k</sup>* <sup>−</sup> 2.38 <sup>×</sup> <sup>10</sup>−4*T<sup>k</sup>* 2 .

#### 2.3.2. Multifactor Attribution Analysis

Multifactor attribution analysis decomposed the impact contribution according to the fixing-changing method and was calculated as follows [27,42].

$$
\Delta X\_j = \frac{1}{2^{n-1}} \sum\_{i=1}^{2^n} a\_{i,j} \times S\_i \quad j = 1, \dots, n \tag{7}
$$

$$A = \sum\_{j=1}^{n} \Delta X\_j = \frac{1}{2^{n-1}} \sum\_{i=1}^{2^n} \mathcal{S}\_n^i \times \mathcal{S}\_i \tag{8}$$

$$\beta\_n^i = \sum\_{j=1}^n a\_{i,j} \tag{9}$$

where ∆*X<sup>j</sup>* is the influence contribution of the *j*th factor, *αi*,*<sup>j</sup>* is the weight coefficient of the *j*th factor corresponding to the *i*th scenario, *S<sup>i</sup>* is the simulated result corresponding to scenario *i*, *n* is the number of factors considered, *A* is the sum of the contributions of all factors, and *β i n* is the sum of the weight coefficients of all factors in the *i*th scenario under the premise of considering n influencing factors.

The contribution rate of each factor to the change in the water cycle elements was calculated with the following equation:

$$\eta\_i = \frac{\Delta X\_i}{\sum\_{j=1}^{n} \Delta X\_j} \quad i = 1, \dots, n \tag{10}$$

#### **3. Results and Discussion**

#### *3.1. Model Calibration and Validation*

The WEP-N model in the SRB was constructed according to a previously published method [28]. The SRB was partitioned into 9544 sub-basins and 29,488 contour zones based on DEM data and river observations. WEP-COR, the predecessor of the WEP-N model, was verified using the stratified soil temperature and liquid water content of the Qianguo Irrigation District and the monthly mean discharge between 1956 and 2000 at the Jiangqiao,

Fuyu, Jiamusi Stations in the main streams of the Songhua River [24]. The WEP-N model was validated using the stratified soil temperature and liquid water content, the daily discharge based on the hydrothermal coupling experiment, and the river flow monitoring experiment during two freeze-thaw periods (2017–2018 and 2018–2019) in the Heidingzi River Basin [28]. The WEP-N model was calibrated and validated using the discharge measured monthly between 1956 and 2018 at the Jiangqiao, Fuyu, and Jiamusi Stations in the main streams of the Songhua River. Data from the Jiangqiao, Fuyu, and Jiamusi Stations were divided into two parts. The data for the period 1956–1990 were used for calibration, and those for 1991–2018 were used for validation. The results of the monthly mean discharge simulations at the Jiangqiao, Fuyu, and Jiamusi Stations are shown in Figure 2. In general, the WEP-N model performed satisfactorily for the SRB and achieved efficiency coefficients of NSE > 0.75 and RE < 5 % for the validation period (Table 1). The simulated flow was suitable for application in the subsequent analyses. *Water* **2022**, *14*, 3170 7 of 15

**Figure 2.** Validation of the WEP-N model at the Jiangqiao, Fuyu, and Jiamusi Stations. dall test method and abrupt changes were detected by the Pettitt test method (Figure 3). **Figure 2.** Validation of the WEP-N model at the Jiangqiao, Fuyu, and Jiamusi Stations.

Jiangqiao 0.80 0.77 4.98 4.85 Fuyu 0.86 0.73 4.26 −0.41 Jiamusi 0.81 0.76 4.51 0.83

The trends in measured annual runoff in the SRB were analyzed by the Mann–Ken-

*3.2. Influence of Climate Change and Water Use on the Annual Runoff Variation in the SRB*

**Table 1.** Validation of the WEP-N model at the Jiangqiao, Fuyu, and Jiamusi Stations.


**Table 1.** Validation of the WEP-N model at the Jiangqiao, Fuyu, and Jiamusi Stations.

#### *3.2. Influence of Climate Change and Water Use on the Annual Runoff Variation in the SRB*

The trends in measured annual runoff in the SRB were analyzed by the Mann–Kendall test method and abrupt changes were detected by the Pettitt test method (Figure 3). The measured annual runoff decreased by 20.3 billion in 63 years A significant abrupt change in the measured annual runoff appeared ca. 1998 (*p* < 0.01). To analyze the impact of climate change and water use on the water cycle, the data series were segregated into two periods using the abrupt change in measured annual runoff ca. 1998 as the dividing line. The data for the period 1956–1998 served as the base period, while those for 1999–2018 served as the change period. The measured annual runoff decreased by 20.3 billion in 63 years A significant abrupt change in the measured annual runoff appeared ca. 1998 (*p* < 0.01). To analyze the impact of climate change and water use on the water cycle, the data series were segregated into two periods using the abrupt change in measured annual runoff ca. 1998 as the dividing line. The data for the period 1956–1998 served as the base period, while those for 1999– 2018 served as the change period.

**Figure 3.** Trend and mutation analyses of measured annual runoff in the SRB. **Figure 3.** Trend and mutation analyses of measured annual runoff in the SRB.

Climate change and water use were the main factors contributing to runoff reduction in the SRB. The influences of climate change and water use on the annual runoff variation were then analyzed. A base scenario (BS) was created to represent the pre-1998 configuration. The meteorological and water use data for the base period were replaced with those of the change period. The other inputs remained unchanged, and four different sce-Climate change and water use were the main factors contributing to runoff reduction in the SRB. The influences of climate change and water use on the annual runoff variation were then analyzed. A base scenario (BS) was created to represent the pre-1998 configuration. The meteorological and water use data for the base period were replaced with those of the change period. The other inputs remained unchanged, and four different scenarios were modeled (Table 2).

BSW Change base period water use data to that of the change period BSM Change base period meteorological data to that of the change period

BSWM Change base period water use and meteorological data to those of the change

Table 3 shows relative changes in temperature, precipitation, and water use. The annual runoff decreased in response to climate change and water use. Compared with the

period

**Scenario Description** BS Base scenario

narios were modeled (Table 2).


**Table 2.** Settings of scenarios for the multifactor attribution analysis.

Table 3 shows relative changes in temperature, precipitation, and water use. The annual runoff decreased in response to climate change and water use. Compared with the BS scenario, the water use was 8.5 billion m<sup>3</sup> higher, the precipitation was 26.0 mm lower, and the temperature was 1.3 ◦C higher in the BSWM scenario.

**Table 3.** Changes in each impact factor.


For the BS and BSWM scenarios, the annual runoff volumes in the SRB were 73.7 billion m<sup>3</sup> and 52.9 billion m<sup>3</sup> , respectively (Table 4). The rate of change in the annual runoff was −28.2% in the BSWM scenario relative to that of the BS scenario. The contribution rates of water use and climate change to the annual runoff reduction in the SRB were 23.0% and 77.0%, respectively. These values were determined by multifactor attribution analysis, revealing that climate change was the dominant factor attenuating annual runoff.

**Table 4.** Contributions of various factors to annual runoff reduction in the SRB.


We analyzed the influence of climate change on the reduction in the annual production flow components. Only climate change sets the scenarios for the comparative analyses based on the BS (Table 5).

**Table 5.** Scenario settings for multifactor attribution analysis.


The annual average temperature increased by 1.3 ◦C while the precipitation decreased by 26.0 mm. Hence, there was a 26.9-mm decrease in annual production flow (Tables 3 and 6).


**Table 6.** Annual production flow influenced by climate change.

We then analyzed the annual production flow components (Table 6). Annual surface, soil, and base flow decreased by 16.7 mm, 0.6 mm, and 9.6 mm, respectively, in the BSTP scenario relative to the BS scenario. The rate of change in the annual surface flow in the BSTP scenario was −12.7% to that of the BS scenario, which was the minimal rate of change in the annual surface, soil, and base flows. Nevertheless, the annual surface flow reduction accounted for 62.1% of the annual production flow reduction. The rate of change in the annual base flow was −36.6%, but the annual base flow reduction accounted for 35.7%. The change in annual soil flow was minimal and accounted for only 2.2%. The observed decrease in total annual production flow in the SRB was caused mainly by the decrease in the annual surface flow. The decrease in annual base flow accounted for a relatively small proportion.
