Unprecedented acceleration of winter discharge of Upper Yenisei River inferred from tree rings

The inflow of freshwater and heat to the Arctic Ocean from rivers is of central importance to global climate variability via impacts on Arctic sea-ice coverage and the ocean's thermohaline circulation (Aagaard and Carmack 1989, Rahmstorf 2002, Fichot et al 2013). Accelerating warming in the Arctic and loss of sea-ice cover unfolding over the last 20 years, known as Polar or Arctic amplification (AA), raise many questions about the feedbacks of terrestrial systems coupled with the Arctic system (Serreze et al 2009, Screen 2014, Prowse et al 2015). Emerging evidence suggests that freshwater originating from the largest Arctic watersheds stimulates the AA trajectory (Yang et al 2014, Prowse et al 2015, Agafonov et al 2016). However, hydrology feedbacks from oceans and sea-ice on the regional basin-scale are only beginning to be recognized due to the short record length and sparse network of hydrological observations (Lammers et al 2007, Landerer et al 2010, Shiklomanov et al 2021).

Eurasia contributes 75% of the total terrestrial runoff to the Arctic Ocean (Shiklomanov et al 2000). Large-scale reanalysis of the hydrological observations and remote sensing assessments show an accelerating rate of discharge for the large Siberian Rivers in recent decades (Syed et al 2007, Yang et al 2007, Shiklomanov et al 2021). An increase of 7% for the period 1936–1999 was first identified by Peterson et al (2002). More recent studies confirm the continuing trend of flow increase for all large rivers of the Eurasian continent and highlight a new historical maximum discharge and soil moisture water equivalent in 2007 (Shiklomanov and Lammers 2009, Tei et al 2013, Holmes et al 2016, Shiklomanov et al 2021). Climate models predict that large flow increases will continue across much of Eurasia through 2090 (Bring et al 2017). Furthermore, the winter flow displays the highest rate of increase, although this can be heavily impacted by flow regulation (Lammers et al 2001, Yang et al 2004, Magritsky et al 2018, Melnikov et al 2019). These trends are not limited to the Eurasian Arctic, as Déry et al (2016) report an 18% increase of river discharge in northern Canada over 1989–2013.

The observed recent changes in the hydrologic regime of the Arctic challenge our understanding of how streamflow responds to climate change. Studies show that intensified precipitation induced by warming climate is a major contributor to the recent rise of annual streamflow, but is not the only factor (Shiklomanov and Lammers 2009, Shiklomanov and lammers 2013). Winter streamflow appears to be particularly sensitive to increasing temperature (Wang et al 2021). Remote sensing data, used to evaluate the linkages between climate and streamflow on the sub-continental scale, reveal large differences in spatial patterns, climate variation, and forcing of permafrost thawing between basins of large Eurasian rivers, such as the Ob, Yenisei and Lena (Troy et al 2012, Wang et al 2021). The spatial distribution of permafrost in the upper reaches of these Siberian rivers varies widely, ranging from isolated (<10% of area underlain by permafrost) for the Ob River and discontinuous (50%–90%) for the Yenisei River to continuous (>90%) for the Lena River (Brown et al 2002). Most notably, the ongoing thickening of the active layer over discontinuous permafrost is contributing most to discharge anomalies (Landerer et al 2010, Makarieva et al 2019). We assume that the spatial and temporal dynamics of permafrost degradation (especially in the discontinuous zone) stimulate the permafrost driver of discharge response for climate change. Extending the record of seasonal discharge helps to gain insights to the longer range of such responses, especially in the past when climate was much cooler. A longer record would better quantify the hydrologic responses across the Arctic watersheds not only to climate change but to direct anthropogenic impacts on the terrestrial runoff as well.

The Yenisei River is the largest contributor of freshwater and energy fluxes of all Arctic rivers (Lammers et al 2007, Shiklomanov and Lammers 2013). The Yenisei River, accounting for 1.5% of global runoff, is 3487 km long, drains a 2580 000 km² area from Mongolia across Siberia, and empties into the Kara Sea. Tree rings can yield a long-term perspective on Yenisei River hydrology and the long-term natural variability of freshwater flux to the Kara Sea. Tree-ring data have long been useful to hydrological studies in arid regions (Schulman 1945). Yet, only recently has this potential been tested in cold regions (MacDonald et al 2007, Agafonov et al 2016, Meko et al 2020) where in some settings trees respond to air temperature variation specifically driven by water levels rather than to precipitation reflecting soil moisture stress. In this study we aim (a) to reconstruct seasonal discharge of the Yenisei River at the upper reaches situated in the discontinuous permafrost, and (b) to examine the long-term variability of the streamflow in response to climate change particularly during the winter season, when discharge trends have been pronounced.

The climate-streamflow relationship is examined with the Water Balance/Water Transport Model (WBM) (Grogan 2016, WSAG 2016). WBM has recently implemented a groundwater (GW) modeling algorithm (MODFLOW), which presently is the state-of-the-art GW modeling framework developed by the USGS (Hughes et al 2017). This model addition is a comprehensive hydrological system that tracks the sources and fates of water as it moves through the hydrological system including glacier melt water, snow melt, rainwater, and quantifying their relative distributions along the river systems and at any location above and below the surface (Grogan et al 2016). WBM has been applied to address a variety of hydrologic questions including global and regional water resource management and water services (Vörösmarty et al 2000, 2005), water availability and ecosystem stresses (Vörösmarty et al 2010, Shiklomanov et al 2016), permafrost impacts on hydrological conditions (Rawlins et al 2003, 2013), land use and hydrologic vulnerability (Douglas et al 2007), water temperature (Stewart et al 2013), and hydro-biochemistry (Wollheim et al 2008, 2015).

2.1. Yenisei river upper reaches

2.2. Discharge data

The Yenisei River at Kyzyl gauge represents a predominantly pristine streamflow regime without any significant regulation and water withdrawal across the watershed. The Kyzyl gauge (code 9002) is located at the confluence of two headstreams called Bolshoy Yenisei and Maliy Yenisei (or Bi-Khem and Ka-Khem in the Tuvan language). The Bolshoy Yenisei rises in the Tannu-Ola Mountain range and the Maliy Yenisei in the Darhat rift valley, Mongolia. The Kyzyl gauge is at 51.72° N, 94.40° E and 615 m asl and has a drainage area of 115 000 km². Seasonal dynamics of streamflow are shown in the hydrograph of figure 1, which highlights the low flow during November to April and the high flow from June to October, with the maximum discharge occurring usually in late June due to snowmelt runoff. The average annual discharge is 1025 m³ s⁻¹. The spring flow doubles that amount (2410 m³ s⁻¹) and in winter the flow is only about one-third of the annual average (285 m³ s⁻¹).

Zoom In Zoom Out Reset image size

The monthly observational discharge data are downloaded from the Regional Arctic Hydrographic Network (R-ArcticNET v.4.0 www.r-arcticnet.sr.unh.edu, Lammers et al 2016). Most recent data are obtained from the State Hydrological Institute, Saint Petersburg, Russia, and have been screened for quality control via calculation of gridded runoff fields across the Arctic region (Lammers et al 2001). The Kyzyl monthly observations were updated for the 1927–2019 period, covering 92 years. Missing data are filled with linear interpolation from nearby gauges of the Yenisei River basin.

2.3. Tree-ring data

2.4. Discharge reconstruction model

Climatic signals in tree-ring chronologies are estimated with the program Seascorr (Meko et al 2011), which examines Pearson correlations and partial correlations between tree-ring index series and monthly data of mean temperature and precipitation aggregated over variable-length seasons (one, three and six months). The correlation analysis uses climate data from weather stations Kyzyl (RSMID#3696, 51.72° N and 94.5° E, 626 m asl) and Abakan (RSMID#29862, 53.77°N and 91.32°E, 254 m asl) for the interval 1925–2018, as well as CRU TS v. 4.04 0.5 degree gridded data for the interval 1901–2019 (https://climexp.knmi.nl, Harris et al 2020).

The statistical model selected for reconstruction was stepwise linear regression of discharge on chronologies (e.g. Woodhouse et al 2006). The model is problematic with just 70 observations (1927–1996) for calibration and a total of 108 potential predictors (36 × 3 = 108) if even a simple lagged model with series at lags t − 1and t + 1 are allowed as predictors of discharge in year t. Protocol for statistical screening of the predictors and model verification is shown in supplementary materials (SM2). For each season, discharge between 1937 and 1996 was regressed stepwise on the 14 selected tree-ring series using p-enter = 0.25 and p-remove = 0.30 and a cross-validation stopping rule. Stepwise regression yielded the reconstruction models for seasonal and annual discharge. The first reconstruction model is

$$\begin{equation}y = a + {b_1}{x_1} + \ldots {b_5}{x_5}\end{equation} \tag{ Model 1 }$$

where $y$ is the reconstructed Nov–Apr discharge in year t; ${x_1}, \cdots ,{x_5}$ are tree-ring indices of mong33, russ258, and CHG in in year t, and of HON and russ230 in year $t + 1$; $a$ is the estimated regression constant; and ${b_1}, \ldots ,{b_5}$ are the estimated regression coefficients. The second model is

$$\begin{equation}y = a + {b_1}{x_1} + \ldots {b_4}{x_4}\end{equation} \tag{ Model 2 }$$

where y is the reconstructed Oct–Sept discharge in year t;${x_1}, \cdots ,{x_4}$ are tree-ring indices of russ227 and CHG in year t, of russ249 in year $t - 1,$ and of SHA in year $t + 1;$ and the regression constant and coefficients are defined as in Model 1. While the stepwise modeling on the 14 potential predictors was restricted to the common period 1927–1996 of those chronologies, the predictors in the final models had somewhat longer common periods. Accordingly, models were re-calibrated using years 1927–1997 for Nov–Apr and 1927–2000 for Oct–Sep to arrive at models applied for long-term reconstruction.

2.5. Water balance modeling

3.1. Hydrological signal of tree rings

3.2. Seasonal discharge reconstruction

Generating a reconstruction from the selected predictor pool was also successful. Both models demonstrate skill in prediction since the reduction of error (RE) statistic is positive in both cases (table 2). There was no indication the model residuals violated any assumptions, as indicated by a non-significant DW statistic (table 2), a normal-looking histogram, and a featureless scatterplot of residuals on fitted values (figure S2). Moreover, a straight line fit of residuals against time indicated no significant trend in residuals. Supplementary materials include the equations of the two reconstructed models: Nov–Apr low water season and annual discharge of Oct–Sept (equations S1 and S2).

Table 2.  Reconstruction statistics for Yenisei River streamflow at the Kyzyl gauge. Calibration interval: 1927–1997 (71 years) in Model #1 and 1927–2000 (74 years) in Model #2.

Predictant	Predictors	R	R2	adjR2	DW	F	RE	RMSEcv
Model #1 Nov–Apr	mong33t russ258t CHGt HONt −1 russ230t +1	0.72	0.52	0.48	1.94 (p = 0.54)	14.14 df = (5, 65) p < 1E-8	0.44	41.2832
Model #2 Oct–Sept	russ227t CHGt russ249t −1 SHAt +1	0.51	0.26	0.22	1.85 (p = 0.40)	6.08 df = (4, 69) p < 0.001	0.11	126.3911

^a R—coefficient of correlation; R²—coefficient of determination; adjR²—coefficient of determination adjusted for number of predictors in model; RE—cross-validation reduction of error; DW—Durbin-Watson statistic and its p-value (p > 0.05 indicates no lag-1 autocorrelation of residuals); F-overall F-statistic for testing significance of regression model, with degrees of freedom and p-value (p < 0.05 indicates significant model); RMSEcv—root-mean-square error of cross-validation, a measure of reconstruction uncertainty.

Stepwise regression of the selected tree-ring chronologies with Yenisei discharge at the Kyzyl gauge finds two seasonal groupings: (a) a 6 month period beginning in November of the preceding growth year and ending in the following April with the strongest signal (R² = 0.52, F = 14.14, p < 0.00001) (table 2), and (b) a 12 month window beginning in October of the preceding growth year with a weaker signal (R² = 0.26, F = 6.08, p < 0.001). Although the signal of the second seasonal grouping was weaker, it is still comparable with signals obtained for the Ob River basin from the Eurasian pan-Arctic (Agafonov et al 2016, Meko et al 2020). Both signals are sufficiently strong to justify the reconstruction effort.

The developed Nov–Apr tree-ring reconstruction tracks the year-to-year variations of baseflow discharge and water-year discharge since 1784. The decadal and multidecadal variability of the two reconstructed flow series is similar and the reconstructions share common features until 1960, after which the agreement falls apart. In particular, the rate of flow in the winter has accelerated considerably over the last 25 years. The winter flow had a prolonged 71 year period (1920–1990) below the historical mean, then suddenly surged up (figure 2, Nov–Apr). The annual flow rises in 1960–1970 then gradually declines over the next decade, stays elevated during 1995–2015 and finally declines again until the present (figure 2, Oct–Sept).

Zoom In Zoom Out Reset image size

Discharge modeled with tree-rings is double the length of the instrumental record at the Upper Yenisei and found (a) the unprecedented upsurge of baseflow from 1995 and (2) a significant disagreement in the rate and trend of fluctuations between the winter and annual flow during 1995–2019.

3.3. Water balance modeling of discharge

The streamflow–climate relationships for two seasonal windows were evaluated using long-term simulations with the water-balance model (WBM at the Upper Yenisei basin for the interval 1905–2019, 114 years). WBM results show high sensitivity of discharge to precipitation and temperature that tracks well the decadal and multidecadal variability of the flow while systematically underestimating the annual variance (figure 3 and table S2). This could be caused by underestimation of precipitation in our input climate dataset especial in the early observations (figure S3), which may poorly represent the high-elevation precipitation pattern. The number of meteorological stations is sparse, and precipitation varies extensively in space in the Sayan Mountains. Moreover, WBM probably does not consider some environmental factors and processes that may regularly contribute to the discharge variance at annual and seasonal scales. For example, due to insufficient climate and land cover data over such a long period, we had to use simplified evapotranspiration calculations and ignore glacier melt and permafrost thaw. This neglects the possible role of permafrost degradation in the recent increase of the flow.

Zoom In Zoom Out Reset image size

Another notable feature of the WBM calculated series is a significant increase of the estimate error in the winter discharge after 1960 (figure 3, Oct–Sept). The climate drivers decreased their rate of change during 1960–2000 while the baseflow continued accelerating higher. It appears the winter flow decoupled from the decadal trends in climate variability during that time. This large discrepancy in the WBM simulated and observed discharge suggests that the feedbacks of the main forcing factors input to the WBM altered over the period 1960–2000. However, the climate–winter flow relationship strengthened after 2000 and the simulated variance closed the gap with the instrumental flow.

Tree-ring reconstruction of seasonal discharge for the Upper Yenisei extends the instrumental record of winter flows at Kyzyl from 1927 back to 1784. This is the first and only reconstructed winter flow in the Siberian region, bringing a new perspective to streamflow dynamics. The reconstruction permits analysis of both the long-term (decadal and multidecadal) and annual variability of discharge. The long-term variability of winter flow differs significantly from the annual flow. The reconstruction finds a remarkable 80% upsurge of the winter discharge over the last 25 years that is unprecedented in the last 235 years, while the annual discharge increases by only 7%. Our reconstructed annual flow corresponds very well to an independently produced tree-ring reconstruction of Selenga River discharge for the interval 1708–1998 (Andreev et al 2016), which also derives part of its flows from precipitation in the mountainous upper headwaters of Altai-Sayan Mountains as well as from arid parts of the Yenisei Basin southeast of Lake Baikal in Mongolia (figure S4).

Many studies generally attribute the recent change of flow rate across the Eurasian pan Arctic to warming air temperature, seasonal changes of precipitation, or snowpack dynamics (Yang et al 2003, 2007, McClelland et al 2004, Wu et al 2005, Zhang et al 2018). However, more recent assessments conclude that the primary driver of the flow change is enhancement of regional surface–groundwater interactions due to degradation of the permafrost (Evans et al 2020). The linear regression slope of air temperature fields (CRU T.4, 1905–2019) is positive overall at the upper reaches of the Yenisei since 1960s (figure S5). Figure S6 shows a positive linear regression slope in active layer depth data for the interval 2000–2018. CRU reanalysis addressed the studied watersheds, finding no extraordinary changes in seasonal pattern of summer temperature after 1960 and 1995 (figure S7). Tree-ring reconstructions of summer temperature and precipitation across the Altai-Sayan Mountains suggest that the rate of precipitation and temperature change in the late 20th century is not unusual in the context of the last 300 years. The late 20th century was relatively cool and wet and the most recent decades warm and dry (Myglan et al 2012, Kostyakova et al 2018, Fedotov et al 2019, Oyunmunkh et al 2019). Strong correspondence of river discharge to seasonal snow cover changes does not fully explain the winter flow rate (Yang et al 2007, Troy et al 2012). Studies covering various time periods show snow cover and snow depth decreased in the studied region between 1936–2007 (Ye et al 1998, Bulygina et al 2009). VIC-simulated seasonal runoff for the interval 1936–1999 indicates opposite trends in snowpack variability between the north and south Yenisei River basin with significant reduction of snowpack in the south (Troy et al 2012). It is evident that neither temperature, seasonal precipitation changes, nor snowpack dynamics are the primary drivers of the documented streamflow change.

It is worth mentioning the role increased fire frequency plays on the degradation of permafrost. Carbon emission modeling finds Siberia to be an increasingly significant region where frequent boreal fires and intense large area burns impact climate and permafrost (Conard et al 2002, Soja et al 2004a, Tchebakova et al 2009, Kukavskaya et al 2016, Gibson et al 2018, Biskaborn et al 2019). Smoke plumes from south Siberian fires (including the Sayan Mountains, especially near Lake Baikal) have been noticeable on the global scale since the late 1990s (Warneke et al 2009). Soja et al (2004b) report on the post-fire impact on the seasonal frozen layer and soil temperature dynamics. After forest fire, soil temperature can increase by as much as 2 °C–6 °C for up to 15 years after the fire, which, as has been observed since the mid-20th century, also increases the depth of permafrost thaw (Kershaw et al 1975, Viereck and Schandelmeier 1980, Van Cleve and Viereck 1983, Furyaev 1996, O'Neill et al 2003, Kirdyanov et al 2020). In connection with forest fires, ground-based data estimates the top of permafrost is subsiding at a rate of 0.52 cm yr⁻¹ in the Yenisei River Basin (Knorre et al 2019). Overall, the increased fire frequency in southern Siberia likely impacts the aquifer level in the discontinuous permafrost and the water-storing capacity of the ground, which enhances ground water discharge to the flow.

Climatologically, the frequent fires over arid south Siberia are linked to a local high-pressure system controlled by the Arctic Oscillation (Kim et al 2020). Anomalous warmth in late winter (0.8 °C–1.2 °C in Feb–March) enhances evaporation, causing earlier ground surface exposure and drier ground in spring that promotes the spread of fire. AA, due to the weakening of midlatitude jets, is linked to the weather patterns that increase the probability of extreme weather in the study region (Cohen et al 2014, Screen 2014). Interestingly, increased flows is both a consequence of AA, due to the same process described above, and a contributor to AA as larger winter flows bring significantly more fresh water to the Kara Sea.

We hypothesize that the recent upsurge in winter flow resulted from melting permafrost, which itself is a consequence of a warming climate. The lack of long-term studies on permafrost degradation at the upper reaches of the Yenisei River complicates our understanding of permafrost's impact on the hydrological regime (Walvoord and Kurylyk 2016). Still, multiple lines of evidence suggest permafrost plays a key role in this upsurge of winter flows. Biskaborn et al (2019) report accelerating permafrost warming in the circumpolar Arctic in the last decade of up to 0.20 ± 0.10 °C in the discontinuous zone. Large-scale satellite gravitational measurements (GRACE) of the Eurasian pan-Arctic water budget attribute the thickening of the active layer to melting ground ice in discontinuous permafrost, directly contributing to the changes of river discharge since 2003 (Landerer et al 2010). Park et al (2016) showed long-term regional trends (1980–2009) of 0.03 m −0.06 m decade ⁻¹ in thickness of the active layer over the permafrost extent, indicating widespread permafrost degradation. Measurements of stable isotope composition of Yenisei water in watersheds with degraded permafrost show that winter flow decouples from the precipitation isotopic signature in the 21st century (Streletskiy et al 2015). Active-layer depth over Siberia steadily increased over 1956–1990 (Frauenfeld et al 2004, Zhang et al 2005) and continued increasing into the 2000s (Luo et al 2016). Finally, the most recent remote sensing data report a rate of 4.4 mm per year over 1980–2016 in thickening of the active layer across Siberia, which increased the storage capacity of discontinuous permafrost and the contribution of ground water to river discharge during the winter months (Wang et al 2021).

These changes have implications for the vegetation that grows there. The boundary between the active layer and permafrost plays an important role as source water storage (Sugimoto et al 2003). Late-summer-fall precipitation infiltrates into the seasonally frozen layer and mixes with ice-thaw water. The impact of increased evaporation of surface water on Eastern Siberian streamflow is not significant (Sugimoto and Maximov 2012), meaning that ongoing thickening of the active layer due to warming spring temperatures implies thermal distress during growth onset.

A warm spring prompts rapid snowmelt, cooler soil temperature and excess moisture, delaying the onset of growth (Agafonov et al 2016, Kannenberg et al 2019). Furthermore, the positive correlation of tree growth with previous summer and fall precipitation suggests that the growth rate of the current year is enhanced by soil-water recharge during the previous fall, the lack of which can limit accumulation of non-structural carbohydrates that provide reserves necessary for spring growth (Choat et al 2018).

Tree-ring proxies of summer temperature and precipitation define the linkage between climate and Yenisei streamflow. Regression of tree-ring width series with discharge is useful for reconstructing seasonal streamflow changes over centuries in arctic watersheds. Tree ring widths account for 52% to 27% of the explained variance of seasonal and annual flow of the Yenisei River. Applying tree rings from high and low elevations allows modeling of at least two, and possibly more, water seasons from a single watershed. Further efforts should be focused toward increasing the sample depth of the tree-ring chronologies to extend the interval of reconstructed flow dynamics back in time. Reconstruction of the Upper Yenisei flow found an unprecedented increase in the winter flow rate between 1995 and 2019, which is not seen in the record going back to 1784. This rate is nearly 80% above the instrumental and reconstructed average. Water balance modeling indicates the weakening of climate and streamflow linkages soon after 1960. Many studies link AA with the intensified degradation of discontinuous permafrost and boreal forest fires, possibly driving the change of winter flow rate in recent decades.

The study is part of the Tree-Ring Integrated System for Hydrology (TRISH project) that is supported by funding from U.S. NSF Polar Office program #1917503 and #1917515. We also acknowledge support from the Russian Science Foundation (RSF) #19-14-00028 to V Myglan. We thank Professor Orlan Oidupaa (Tuva State University) for the generous assistance in the field, and Alec C Gagliano and Misha Seroukhov (students of University of Arizona) for the lab assistance.

The data that support the modeling and findings of this study are available from the corresponding author or the University of Arizona Library (https://doi.org/10.25422/azu.data.17086889).

All data that support the findings of this study are included within the article (and any supplementary files).

I P P conceptualized the study, developed the tree-ring datasets, modeled and analyzed the Yenisei discharge, and wrote the manuscript. D M assisted in computer programming, modeling and in editing the draft. A S H provided hydro-climatic instrumental data and simulated W B M series. R T discussed and edited the draft, helped with figures drawings. V M, V V B and A V T collected the tree ring samples in Tuva and crossdated the tree-ring materials from S T A and C H G sites.