The present study was conducted in the LYL basin (53°23′53″–53°27′30″N, 122°14′27″–122°21′02″E), which is located at the southern margin of the high-latitude permafrost region in the northernmost China and is a primary tributary of Heilongjiang watershed (Figure 1). The area of the study basin is 21.9 km², and the elevation ranges from 302 to 696 m above sea level, with an average slope of 13.8°. The research zone is characterized by a cold -temperature continental monsoon climate, with a duration of winter for more than 7 months. The annual average temperature is −4°C. The average annual precipitation is 500 mm, of which about 3/4 occurs from June to September. The soil type in the study watershed is dominated by brown coniferous forest soil with the average thickness ranging from 15 cm in the hillslope to 50 cm in the flat valley. The study area also features a permafrost layer. The main tree species are composed of 30% Larix gmelinii (Rupr.) Kuzen., 30% Betula platyphylla Sukaczev, 20% Populus davidiana Dode and 20% Pinus sylvestris var. mongholica Litv.

It is important that higher sampling frequency can provide a more accuracy estimations of young water fraction. For example, compared to weekly sampling, daily sampling resulting in a higher Nash-Sutcliffe Efficiency and more pronounced short term dynamics in the simulation result (Stockinger et al., 2016). In this study, the study period is from June 8^th to October 19^th in the year of 2019. The rainfall was manually sampled, and when multiple rainfall events occur on a day, the δ value was calculated as the weighted average of the rainfall. The runoff was manually sampled at a fixed daily sampling time of 18:00. All samples were collected in 50 mL polyethylene bottles and immediately sealed with Parafilm. Return to the laboratory for filtration and sealing, keep frozen to prevent evaporation, and tested within a few days. Total of 21 rainfall samples and 133 runoff samples were collected during the 133 days study period.

The instantaneous water level height and conductivity data of runoff were monitored by the water level gauge (Onset HOBO U20-001) and conductivity gauge (Onset HOBOU24-01) placed on the composite weir. In order to obtain continuous flow data, a water level flow relationship curve was established based on the known measured water level and the composite weir cross-section flow. The dynamic changes of soil temperature (Ts, °C) and soil moisture (Ms, %, volumetric soil moisture content) observed at fixed points in the watershed using the soil moisture and temperature sensor (Campbell Scientific CS650) at depths of 5, 10, 20, and 40 cm, with recording intervals of 30 min. Rainfall (R, mm) and air temperature (Ta, °C) data were provided by meteorological station immediately adjacent to the watershed (Figure 1). The dynamics of Ta, Ts, R, and Ms in the study period can be found in Figure 2.

All isotope data were analyzed at Shenyang Institute of Applied Ecology using the liquid water isotope analyzer (LGR-DLT-100, United States). The stable isotope composition of ¹⁸O and D was expressed in δ (‰ or per mil), with its definition based on the Vienna Standard Mean Seawater (VSMOW) measurement of δD and δ¹⁸O (where δ(‰) =(R_sample/R_standard-1)×1000). The analytical accuracy of δ¹⁸O and δD are ±0.10‰ and ±0.30‰, respectively.

The introduction of water isotope (¹⁸O or D) as a tool to divide the components of event water and pre-event water based on their source time in hydrograph separation is because isotope-based hydrograph separations (IHS) rely on the stable water components that can be objectively measured. In the mass balance method can be used to separate the runoff components if the initial isotope values of the two water sources are significantly different. The calculation procedure can be defined as follow: Q t = Q p + Q e (1) C t Q t = C p Q p + C e Q e (2) F p = C t ‐ C e C p ‐ C e (3a) where Q _t represents the runoff depth (mm/day), Q _p represents the pre-event water (mm/day), Q _e represents the event water (mm/day), C _t , C _p and C _e represent the δ values of runoff depth at the timing of t, pre-event water and event water (‰), and F _p represents the fraction of pre-event water in the stream.

Theoretically, the pre-event water is the water stored in the basin before the event occurs. The pre-event water flows through the soil and rock formations and then flows into the river channel together with the event water to form river runoff. Thus, the ion concentration (conductivity) of the pre-event water is typically greater than that of event water. The conductivity mass balance method (CMB) also relies on the mass balance method to identify the proportion of different water sources. Where Q _t is the runoff depth (mm/day), Q _p is the component of pre-event water (mm/day), Q _e is the component of event water (mm/day), CD _t , CD _p and CD _e are the conductivity values of stream water at the timing of t, pre-event water and event water (μS/cm), and FD _p is the fraction of pre-event water in runoff using CMB method with the equation as follow: F D p = C D t ‐ C D e C D p ‐ C D e (3b)

In this study, we define pre-event water (usually called old water) as the water stored in the catchment before the rainfall event, while the event water (usually called new water) is the water from the current precipitation event.

Young water fraction (F_yw) represents the proportion of water that is less than a certain threshold age. The F_yw of the study catchment was calculated using the method developed by Kirchner. (2016a). First of all, sine wave fitting was performed on the time series of δ¹⁸O (‰) values for rainfall and streamflow, and determine the sine and cosine coefficients of rainfall and streamflow through multiple linear regression: C r t = a r ⁡ cos 2 π ft + b r ⁡ sin 2 π ft + k r (4) C s t = a s ⁡ cos 2 π ft + b s ⁡ sin 2 π ft + k s (5)

Where C _r (t) and C _s (t) are the actual measured δ¹⁸O (‰) values of rainfall and runoff, respectively, and a and b are the coefficients used to determine the amplitude and phase shift of the seasonal δ¹⁸O cycle, f is the frequency of the period (f=1 year^-1 for a seasonal period), t is the time (in days) after the start of the sampling period, k _r and k _s are constants for sine wave fitting: A r = a r 2 + b r 2 , A s = a s 2 + b s 2 (6) φ r = tan ‐ 1 a r b r , φ s = tan ‐ 1 a s b s (7)

Where A _r and A _s represent the amplitude of rainfall and runoff, respectively, φ _r and φ _s represent the phase shift of rainfall and runoff, respectively.

The phase shift difference can be expressed by the shape parameter α in the gamma distribution function Γ(α,β) as: φ s ‐ φ r = α tan ‐ 1 A s / A r ‐ 2 / α ‐ 1 (8)

The scale parameter β in Γ(α, β) can be expressed as: β = 1 2 π f A s / A r ‐ 2 / α ‐ 1 (9)

The threshold age τ _yw of young water fraction can be obtained by using the regression equation proposed by Kirchner (2016a): τ yw / T ≈ 0.0949 + 0.1065 α ‐ 0.0126 α 2 R 2 = 0.9998 (10)

F_yw can be obtained by lower incomplete gamma function Γ (τ, α, β). F yw = Γ τ yw , α , β = ∫ τ = 0 τ yw τ α ‐ 1 β α Γ α e ‐ τ β d τ (11)

MTT: mean transit time, the average time for water parcels to enter the catchment in the form of precipitation until it leaves the catchment. It can be expressed as: MTT = 1 ‐ C 2 π C (12)

Where C is the attenuation coefficient, C=A _s /A _r .

In the current study, Pearson correlation analysis was applied to determine the relationship between the δ¹⁸O of the runoff and Ta, rainfall, Ts and Ms. Subsequently, in order to eliminate the effect of multicollinearity on the equation, the arithmetic mean values of Ms and Ts from different layers (depths of 5, 10, 20, and 40 cm) in this article were selected. Finally, a multiple linear regression including δD and δ¹⁸O of the runoff was performed to identify the main influencing variables. The backward trajectory calculation was obtained by Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) web version combined with ArcGIS. The fitting of the sine wave yielded the F_yw of the rainfall and runoff by using Matlab.

The δ¹⁸O value of rainfall in the LYL watershed varied between −25.0‰ and −4.5‰, with an average value of −13.1‰, and the δD value varied between −188.2‰ and −40.9‰, with an average value of −100.3‰. Furthermore, the δ¹⁸O value of the river water varies between −16.8‰ and −14.1‰, with an average value of −16.2‰, and the δD value varies between −124.6‰ and −112.5‰, with an average value of −119.0‰.

Figure 3 shows the relationship between the δD and δ¹⁸O for rainfall and stream water in the LYL watershed, along with a comparison to the Global Meteoric Water Line (GMWL). Craig (1961) proposed GMWL for global rivers, lakes and rainfall samples, has a relationship of δD=8δ¹⁸O+10‰. The Local Meteoric Water Line(LMWL) (δD = 7.41δ¹⁸O-3.05‰, n = 21, R ²=0.97, p < 0.001) of rainfall in the LYL watershed has a similar slope and a smaller intercept compared to the GMWL. However, there is a significant difference between the LMWL (δD=4.82δ¹⁸O-41.00‰, n=133, R ²=0.61, p < 0.001) of the stream and GMWL, with a gentler slope and a larger intercept. The LMWL slope (4.82) of the runoff is the lowest than for the LMWL of rainfall and GMWL.

Comparing the IHS and CMB (Figure 6), it can be seen that the dynamic of the proportion of new water in runoff components is consistent. Both methods demonstrated that runoff mainly derived from old water (pre-event water). We also noted that CMB underestimated the pre-event water fraction during the rainfall period compared to the isotope method. Additionally, we found that CMB was not able to capture the minor changes of the runoff components when an event with small rainfall amount occurred. For example, from early to mid August, CMB shows that the runoff was composed of old water and does not reflect the proportion of new water. HIS could divide small proportion of new water.

Isotopic hydrograph separation analysis showed that the contribution of new water and old water was 3.6% (Table 2) and 96.4% during the study period, respectively. The rainfall amount and duration of the first peak in summer (Event I) were 127 mm and 9 days, respectively (Figure 7). Event I is composed of 7.5% new water and 92.5% old water. The maximum contribution of new water at the peak is 13.0%. The rainfall amount and duration of Event II (the second summer peak) were 123 mm and 7 days respectively. Event II is composed of 12.5% new water and 87.5% old water. The maximum amount of new water in the runoff component at the peak reaches 24.5%. Furthermore, the runoff calculated by conductivity mass balance method is composed of 5.7% new water and 94.3% old water during the study period. Event I is composed of 15.4% new water and 84.6% old water. The maximum contribution of new water at the peak of the flood is 34.1%. Event II is composed of 20.8% new water and 79.2% old water. New water accounts for the largest proportion of the streamflow, reaching 48.8% at the end of August (Table 2).

It can be seen from Figure 7 indicated that the rainfall in Event I and Event II had similar rainfall amount and duration. Interestingly, we found that the shift in soil moisture from dry to wet states significantly influenced the composition of streamflow (Figure 7). Specifically, the soil moisture content before Event I was 0.172%, while the soil moisture content before Event II was 0.250%. The higher the soil moisture before rainfall lead to the greater the new water fraction. Moreover, the degree of underestimation of the old water fraction using the electrical conductivity method increased with an increase in the soil moisture before rainfall.

The young water fraction (F_yw) of the LYL watershed was estimated using the rainfall and runoff δ¹⁸O values measured by high-frequency daily sampling data. The coefficient of sine and cosine through multiple linear regression was determined by performing sine wave fitting on the δ¹⁸O time series data of rainfall and runoff and calculating the amplitude of the sine wave fitting (Figure 8). The largest deviation of rainfall and runoff from their corresponding annual average isotope content was related to the estimated F_yw (Table 2). The sine wave fitting of rainfall and runoff was found to be statistically significant (p < 0.0001). The F_yw value of the runoff was 5.6% for the LYL watershed. According to the above formula, it can be calculated that, the threshold age is 65 days (Table 3). The uncertainty ranges of F_yw and τ_yw were 1.8%–24.5% and 58–82 days, respectively. The MTT was equal to 3.33 years.

The isotopic composition of rainfall and its variation are the most important basic data for studying river runoff, lake water and groundwater. Due to the isotope fractionation in the process of evaporation and condensation (Craig, 1961), the hydrogen and oxygen isotopic composition of atmospheric precipitation has a linearly related change, that is, the atmospheric waterline. The slope of LMWL is 7.41 less than GMWL, which means that the study area is in the inland with low temperature and far away from the vapor source. The intercept of LMWL is smaller than GMWL, which indicates that the water vapor has experienced multiple condensations during the process of reaching the basin, and the remaining water vapor is formed δD in rainfall is relatively depleted and deviates from equilibrium (Araguás-Araguás et al., 2000; Wang et al., 2019; Keesari et al., 2021). This is similar to the atmospheric waterline relationship obtained in northern China (Zhao et al., 2018). Compared with the LMWL of rainfall and the GMWL of rainfall, the stream water LMWL has the smallest slope and the largest intercept. On the one hand, it shows that the stream water is mainly recharged by the pre-event water (Marc et al., 2001). On the other hand, the stream water is most affected by evaporation and enriched in isotopes (Sanchez-Murillo et al., 2015). Compared with GMWL, a smaller slope represents relatively strong evaporation, which is similar to the results of river research in cold temperature regions of China (Li et al., 2014).

The entire sampling period was divided into monthly intervals from June to October, and the backward trajectory of each rainwater sample was calculated for 5 months (Figure 9). The water sources of the LYL watershed include four regions. In June, it mainly comes from Siberia. In July, in addition to Siberia, there is also part from the Arctic Ocean, and a small part comes from the yellow sea and South China Sea. In August, in addition to Siberia, there is also part from the Pacific Ocean. It mainly comes from Siberia and the Arctic Ocean in September, and Siberia in October. Similar conclusions were reached with previous studies (Li et al., 2012; Ding and Gao, 2015). The changes in time series isotope values fully confirm the variability of water vapor sources and the complexity of evaporation and condensation experienced in the transportation process. Comparing the rainfall δ¹⁸O in June, July, August and September and October, it is found that the δ¹⁸O in September and October is more depleted because of the same tracking 10-day trajectory, which is comes from cold regions, and the longer transport path is more prone to multiple condensation events.

Comparing the results of IHS and CMB hydrograph separation, although IHS sampling is cumbersome and expensive, it is generally considered to be more accurate (Mcdonnell et al., 1990; Laudon and Slaymaker, 1997; Klaus and McDonnell, 2013). The overall trend of the two methods is consistent with the key nodes, so CMB, which can obtain continuous data and has a lower cost, is also a good substitute. To obtain more accurate and detailed results, CMB is not optimal compared to HIS. Because CMB not only underestimates the components of the old water when a rainfall event occurs, but also fails to respond and feedback in time when a small amount of rainfall occurs. It may be because the conductivity of the medium itself is affected by too many external factors and cannot reflect subtle changes. δ¹⁸O is a natural tracer and extremely sensitive, which can reflect extremely subtle changes in the substance itself.

Comparing Event I and Event II, it can be found that although the rainfall and rainfall duration of the two rainfall events are similar, the runoff depth and runoff composition (whether by IHS or CMB method) are significantly different. The soil moisture condition before the rainfall event is the root cause of this difference (Penna et al., 2011; Pan et al., 2018). Before the occurrence of Event I, the soil was in a relatively dry state. At the initial stage of Event I, there was less runoff and more rainfall was given priority to replenish the soil water deficit, which would lead to a lag in the response time of runoff to rainfall (Penna et al., 2011). There were 6 rainfall events in the 9 days before Event II, which increased the area of the saturated zone in the basin and enhanced the water conduction capacity, and the rainfall could be quickly converted into a full runoff (Penna et al., 2011; Penna et al., 2016). So it can be concluded: the greater the soil moisture before rainfall, the greater the new water fraction.

Considering comparing F_yw and MTT in calculating water age at the same time, Hu et al. (2020) explained that lower F_yw and longer MTT means that the response time of water resources management has a considerable lag, which is meaningful for the evaluation of river water flow and water quality. The underlying surface of the karst system is complex, and comparing F_yw and MTT can track and quantify the flow path very well (Rusjan et al., 2019). Our research also selected both the estimation of F_yw and MTT.

This is the first time that the F_yw of a small forest watershed has been estimated in the Daxing’an Mountains, and the results obtained are consistent with the global rainfall and runoff isotope database compiled by Jasechko et al. (2016). Our research concluded that the young water score of 5.6% is lower than the global average of 26%. The first reason may be because the study basin is located in a permafrost area, and the flow of permafrost rivers is mainly controlled by the melting and freezing of the active layer (Wang et al., 2009), that is, when the active layer melts, more old water is quickly supplied to the runoff. Song et al. (2017) also illustrates this point in the permafrost regions of the Qinghai-Tibet Plateau. The second reason may be due to the strong capacity of soil water conservation in small forest watersheds (Wang et al., 2013). The larger the young water component means that it can accelerate the biological cycle process and accelerate the transfer of pollutants, thereby affecting the entire aquatic system (Jasechko et al., 2016).

Some small watershed studies have shown that MTT is generally less than 1 year (Tetzlaff et al., 2010; Muñoz-Villers et al., 2016; Rusjan et al., 2019). The 3.33 years obtained by our research were longer than those of the previous studies, probably because of the environmental attributes of the frozen soil area in the small forest watershed, and the change of the active layer will cause the change of MTT. Because slope, precipitation, snowmelt, and canopy closure are all the major controlling factors that affect the time and size of water flowing through mountainous basins (Wolock et al., 1997; Tetzlaff et al., 2010; Stockinger et al., 2017; Fang et al., 2019).

As a fragile ecosystem, these findings are of great significance to the development, utilization and protection of water resources. Our estimation of young water components and mean transit time based on isotope method is helpful for enhance understanding of the hydrological process, biological cycle, and effective protection of water resources in permafrost regions. As a sensitive permafrost region, the lower F_yw value and the higher MTT value indicate that the active layer stores and provides a lot of water resources for the stream.