Influence of snowmelt on increasing Arctic river discharge: numerical evaluation

Snow is the most important component of the Arctic climatic and hydrological system and is directly vulnerable to climate change. In recent decades, observations have indicated significant decreases in the Arctic snow cover and snowfall rate, whereas water discharge from circumpolar Arctic river basins into the Arctic Ocean has increased. To evaluate the contribution, not well quantified, of snow to the river discharge increase, we conducted sensitivity simulations with surface air temperature and precipitation as climatic treatment variables, combining a land surface model and a distributed discharge model. Variables were treated assuming higher climate variations in the Arctic cold season in 1979–2018. The surface and subsurface runoffs simulated by the land surface model were set as inflows in the discharge model to estimate river discharge. Snowmelt mostly converted to surface runoff, accounting for 73.6% of the anomalous surface runoff increase and inducing the simulated peak discharge in spring and early summer. This relationship was enhanced by the winter precipitation increase. Snow loss induced by higher air temperature contributed to the decrease in the peak and annual discharges, but caused the peak discharge to occur earlier. Additionally, warmer temperature increased the proportion of rainfall in the partitioning of precipitation, causing more subsurface runoff, particularly in autumn and winter. These results provide a first separate evaluation of factors influencing Arctic water discharge, including seasonal hydrographs, and illustrate the influence of climate warming-induced snowfall and rainfall variations on the circumpolar Arctic river discharge.


Introduction
The Arctic is a key region for climate change, where environmental changes are the most drastic and air temperature is increasing three times faster than the global mean (IPCC 2022).However, despite its vulnerability to climate change, snow cover, which is closely linked to climate and hydrology, remains complete in Arctic circumpolar regions during the cold season.Extensive observations have indicated considerable changes in the snow conditions in recent decades.For example, a retreat of the snow cover at the pan-Arctic scale (Mudryk et al. 2021;Robinson 2021) and an earlier onset of melting (Mudryk et al. 2017) have been observed in spring.Moreover, the snow water equivalent (SWE) in winter has decreased, particularly in North America (Pulliainen et al. 2020) where snow cover variations exhibited regional differences, with a decreasing trend in Alaska but an increasing trend in eastern Canada (Hessilt et al. 2023).These changes are connected to modifications of the circumpolar Arctic river hydrograph, because its spring peak-flow regime is strongly related to snowpack melt.In fact, earlier snowmelt resulted in a notably earlier river discharge peak flow in spring (Tananaev et al. 2016;Suzuki et al. 2020) and an increase in the thermal energy inflow at the Page 2 of 15 Park et al. Progress in Earth and Planetary Science (2024) 11:13 soil surface, causing earlier soil thawing (Kim et al. 2012;Park et al. 2015) and further adding snow meltwater to soil water storage.Consequently, increasing soil moisture in the thawing active layer enhanced evapotranspiration and runoff (Park et al. 2020(Park et al. , 2021(Park et al. , 2022)).In circumpolar regions, evapotranspiration largely consumes the water of the upper soil layer, mixed by snow meltwater and rain water, owing to the relatively shallow plant root depth (Park et al. 2021).In summer, evapotranspiration is greater than precipitation (Park et al. 2008), possibly leading to seasonally drier soil (Hiyama et al. 2023;Park et al. 2022).Additionally, owing to the warming climate, evapotranspiration has shown a trend of increase in recent decades (Helbig et al. 2020;Kim et al. 2023).Enhanced evapotranspiration can weaken the link between snow meltwater in soil water storage and summer rainfall to runoff, consequently reducing summer river discharge, as observed in the case of Arctic rivers (Déry et al. 2016;Shiklomanov and Lammers 2009).Linked to relatively low evapotranspiration in autumn, increased autumnal precipitation (Box et al. 2019) contributes to greater soil water storage and thus higher runoff during winter (Hiyama et al. 2023).Additionally, the recent warming has caused delay in the onset of soil freezing (Kim et al. 2012;Park et al. 2016), lengthening the period over which soil water is likely connected to runoff (Evans et al. 2018).During winter, snow cover functions as an insulator for the permafrost soil thermal regime.Anomalously thick snow insulates permafrost not only during the winter season, but until the following summer, thereby contributing to the deepening of the active layer (Park et al. 2013(Park et al. , 2015)).Active layer development enhances the melting of permafrost ground ice (Hiyama et al. 2023), potentially strengthening subsurface runoff (Evans et al. 2020;Han and Menzel 2022;Walvoord and Striegl 2007;Wang et al. 2021).Such complex hydrological processes, strongly affected by snow cover variability, likely cause anomalous river discharge in the Arctic.However, the river discharge effectively measured in North America exhibited an increasing trend in the last decades (Feng et al. 2021), inconsistent with the observed negative SWE trend (Pulliainen et al. 2020).Because discharge is the balanced output of precipitation, evapotranspiration, and soil moisture, its observed trend cannot be easily correlated with snow variability.However, the spring river flow regime, strongly dependent on snowmelt, likely reflects snowmelt anomalies.Indeed, a previous study reported that variations in the maximum daily discharge and cold season precipitation were consistent for most Russian rivers (Shiklomanov et al. 2007).
As mentioned above, the observed snow cover in Arctic regions reflects regionally different trends in recent years (Pulliainen et al. 2020).The heterogeneity is likely linked to regional anomalies in snowmelt-associated river discharge, at least in the period of peak flow in spring and early summer.At the onset of snowmelt, river ice also begins to break, and ice floes combine with runoff and the rapid flood pulsing from surrounding uplands to cause ice jam flooding (Rokaya et al. 2018).This phenomenon has huge impact on biogeochemical and socioeconomic systems on local-regional scales (Prowse and Beltaos 2002).Similarly, the impact of snowmelt on river discharge and on human society has been widely recognized.Previous studies established that snowmelt accounted for approximately 60% of the total river discharge for the Kolyma River Basin (Welp et al. 2005) and Yukon River Basin (Lachniet et al. 2016).However, such assessment needs verification using additional observed tracer data such as stable water isotopes in river discharge.If heterogeneous climate features and the geomorphology of the observed rivers are considered, then the observed data are limited in terms of expansion to other Arctic rivers.Consequently, there are no similar assessments for other Arctic regions to apportion the absolute volume of river discharge inflowing to the Arctic Ocean.These limitations could be partially mitigated through application of numerical modeling.
Additionally, the contribution of snowmelt to river discharge is influenced by climate variations; specifically, it decreases with increasing air temperature.For example, the observed Arctic terrestrial mean air temperature in the Arctic cold season (October-May) has increased by 3.1 °C from 1971 to 2017 (Box et al. 2019), possibly enhancing snow sublimation in winter and evaporation in spring, and consequently reducing the fraction of snowmelt contributing to river discharge (Suzuki et al. 2015).Under climate warming conditions, this decrease could be strengthened by increasing rainfall rates (Bintanja and Andry 2017).Most climate model projections indicate stronger warming in the near future (Tebaldi et al. 2021), which would alter the relationship between snow and river discharge (Berghuijs et al. 2014).
In this study, we asked the following questions: • How much does snowmelt currently contribute to river discharge in the Arctic?• In what regions is river discharge most sensitive to snow variability?• How do climate warming-induced snowpack variations influence river discharge?
To answer these questions, we quantify the sensitivity of the circumpolar Arctic river discharge to snow cover variability in 1979-2018 during the Arctic cold season, using simulations from coupled land surface and distributed discharge models.Model simulations are conducted for several scenarios, depending on two climatic treatment variables (surface air temperature T a and precipitation P).Results for important diagnostic variables are compared against a reference simulation with presentday climate conditions.Thus, anomalous output values indicate the influence of snow cover or, more generally, climate variability, on river discharge.Because this study focuses on the sensitivity of river discharge to climate change-induced snow variability, not on its amplitude, evaluation of model results against observations is not included in this analysis.

Model description
In this study, we used the process-based "coupled hydrological and biogeochemical" (CHANGE) land surface model, based on hydrological, physical, and biogeochemical principles (Park et al. 2011(Park et al. , 2018)), to evaluate the impact of snow and climate on the circumpolar Arctic river discharge.CHANGE simulates the heat, water, and carbon dioxide fluxes in the Arctic terrestrial atmosphere-vegetation-snow-soil system.The snowpack preserves snow water and the rain water resulting from higher temperatures during winter.The snow meltwater in spring is divided into soil infiltration and surface runoff, depending on water saturation and freezing/thawing conditions in the soil surface layer, even during precipitation events.With snow meltwater, the excess water from the surface layers of the soil, frozen during winter, is converted to surface runoff, which is defined as snowmelt in spring.The infiltrated water flows down or up, depending on the water potential gradient between soil layers.Excess water in either the permafrost surface layer or the soil bottom boundary layer generates a subsurface flow.The soil is dried by surface evaporation and plant transpiration and wetted by precipitation.It experiences phase changes in the Arctic circumpolar region, seasonally alternating between freezing and thawing.Freezing/ thawing phase changes in a soil column (depth: 50.4 m) are an important physical process in the model to explicitly represent the dynamics of water and heat fluxes in both permafrost soil and deeper ground.
The CHANGE-simulated surface and subsurface runoffs are set as inflows for the storage-based, distributed river routing "total runoff integrating pathways version 2" (TRIP2) model (Ngo-Duc et al. 2007) to simulate river discharge and storage (Park et al. 2016(Park et al. , 2017)).TRIP2 represents a more realistic travel time without requiring parameter calibration from observational data.It calculates riverine water storage and river discharge, in each grid cell and at every time step, with the global river network TRIP (0.5° spatial resolution; Oki and Sud 1998).Water storage in individual grid cells is calculated prognostically from the inflow-outflow balance; other variables are diagnosed from the water storage term.Surface and subsurface runoffs simulated by CHANGE are directly stored into individual river storage reservoirs in TRIP2, and then water flows along a prescribed channel network.The TRIP2 model also includes a river-ice algorithm that simulates seasonal ice formation, growth, and melt, depending on temperature variability of the river water (Park et al. 2016).River water temperature along the river network is calculated with a one-dimensional advection-diffusion equation, considering inflows of heat and water, both from upstream and from tributaries, and heat exchange at the air-water interface (Park et al. 2017(Park et al. , 2020)).In the calculation of heat exchange with the atmosphere, river storage, and ice volume, the river surface area represents a critical variable that is estimated using two terms: the river channel length within the grid box, calculated geomorphologically based on the TRIP river network data, and the river width, estimated using a geomorphological relationship with mean annual runoff (Arora and Boer 1999).Reservoir regulation, affecting downstream discharge in some river basins, was not considered in our simulations.

Data
Precipitation is an important meteorological variable to simulate river discharge.Arctic regions are generally characterized by strong winds and sparse weather stations, inducing large measurement biases in the observed precipitation.Thus, using a forcing dataset derived from such data increases uncertainty on the simulated river discharge.To minimize this uncertainty, we selected a gridded climatic dataset as input for the CHANGE and TRIP2 simulations: the bias-corrected reconstruction of near-surface meteorological variables from the latest European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalyses (ERA5), obtained with the same methodology as the water, energy and climate change (WATCH) forcing data.This climatic dataset is known as the "WATCH Forcing Data methodology applied to ERA5" (WFDE5, https:// doi.org/https:// doi.org/ 10. 24381/ cds.20d54 e34).The WFDE5 dataset is provided with a daily time step, aggregated from hourly-timescale data, from 1979 to 2018, on a regular latitude-longitude grid with a spatial resolution of 0.5° × 0.5°.The original horizontal resolution of the ERA5 data product (0.25° × 0.25°) has been degraded in WDEF5 to facilitate comparisons and bias correction, then corrected for elevation and for monthly biases derived from the observational dataset of the Climatic Research Unit of the University of East Anglia, England (Cucchi et al. 2020).Finally, using the adjusted WFDE5 data as input for TRIP2, we simulated river discharge and water temperature in each grid cell of the circumpolar Arctic river basins.The circumpolar Arctic river drainage system and river basin boundaries are illustrated in Fig. 1.

Sensitivity simulations
The WFDE5 dataset includes the daily mean, maximum, and minimum values for T a , P, water vapor pressure, specific humidity, wind speed, and solar radiation.With these data, we simulated the water budget in circumpolar Arctic river basins for the 40-year study period.In our simulations, land cover was fixed to its present-day distribution (Lawrence and Chase 2007), while vegetation phenology was simulated prognostically at every time step from the calculated carbon and nitrogen contents, climate conditions, and soil temperature and moisture.Soil thermal and hydraulic properties were calculated with CHANGE using soil texture vertical fractions for sand, silt, and clay, combined with the simulated soil organic matter content, at hourly time steps and with the nominal WFDE5 spatial resolution.
To diagnose the influence of T a and P variations during the Arctic cold season (i.e., September to May) on snow cover and, consequently, on river discharge and soil thermal and hydrological regimes, we designed ten scenarios for our numerical simulations.The simulations comprised a reference simulation, using the original WFDE5 product as forcing data, and three sets of three treated simulations, also forced by WFDE5 data, but with additional anomalous increases in T a , P, or both simultaneously (Table 1).The mean value of P in circumpolar Arctic river basins was taken as 46.5, 29.6, and 31.9 mm in September-November, December-February, and March-May, respectively.Anomalous T a and P increases in the sensitivity simulations were derived from previously published results, briefly described here.T a is a critical variable for soil freezing/thawing and thus temperature and snow cover variation, useful in assessing both warming trends in the current climate (Box et al. 2019) and future warming projections (IPCC 2022).In the Arctic, increases in terrestrial T a of 1.5-2 °C over the previous four decades have been reported (Bekryaev et al. 2017).Increasing Arctic T a also induces P increases.In the central Lena River basin in Siberia, large P increases in 2004-2008, by 30% or more relative to the 1979-2013 average, have been measured (Iijima et al. 2016).In this study, we used these large variation rates as thresholds to define higher climatic changes, then conducted model simulations with T a and P as treatment variables to assess snowpack sensitivity to higher climate warming and the resulting river discharge.Using such climate thresholds in the model simulations allowed us to evaluate the sensitivity of target diagnostic variables to the treated meteorological variables (Park et al. 2015).In the reference simulation, no treatment was applied.In the first set of three simulations (EX TA_SON , EX TA_DJF , and EX TA_MAM ), only T a was treated, with a constant increase of 2 °C in different periods of the Arctic cold season, consistently applied every year throughout the study period.In the second set of simulations (EX PR_SON , EX PR_DJF , and EX PR_MAM ), only P was treated, with an increase of 30% in the same three periods as for the first set.Finally, the third set (EX TA_PR_SON , EX TA_PR_DJF , and EX TA_PR_MAM ) Fig. 1 Study area: the circumpolar Arctic river drainage system.In this region, we conducted model simulations to analyze Arctic river discharge during the study period .Boundaries of the circumpolar Arctic river basins were derived from the sea boundaries of Lammers et al. (2001)  During precipitation events in the cold season, P partitions into snowfall (P s ) and rainfall (P r ) depending on T a .In the CHANGE model, partitioning is achieved relative to a temperature threshold, as defined by Wigmosta and Lettenmaier (1994): where T max is the maximum air temperature at which snow occurs and T min is the minimum temperature at which rain occurs.Because T max and T min were fixed to 0.5 °C and − 1.0 °C, respectively, in our model simulations, partitioning of P into P s and P r largely depended on T a variability.Thus, these equations indicate that the influence of a T a increase on P partitioning should be particularly measurable in spring and autumn, the transitional seasons between cold and warm climate conditions.

Snowfall and rainfall
In the reference simulation, annual mean P s and P r derived from the WFDE5 product (spatial distribution of the annual mean trends, Fig. 2a, b, respectively; interannual trend variability, Fig. 2c), averaged over all grid cells of the circumpolar Arctic river basins (Fig. 1), accounted for 42% and 58%, respectively, of the annual mean P (Fig. 2c).Annual mean P s exhibited a small increasing trend, not statistically significant (0.03 mm year −1 , p < 0.41).Conversely, P r exhibited a strong significant increasing trend (0.87 mm year −1 , p < 0.01).There were also clear differences between P s and P r in the spatial distribution of their annual mean trends.Annual mean P s trends showed values consistently within ± 3 mm year −1 , without large regional differences except in a few small regions with larger decreasing trends (mainly easternmost Russia and southwestern Norway, Fig. 2a).Overall, regions with positive (but not statistically significant) P s trends accounted for 58.2% of the combined circumpolar Arctic river basin area, markedly more than the negative trend regions (41.9%).On the contrary, P r displayed strong regional discrepancies, with strong increasing trends of up to 8 mm year −1 on the Eurasian continent, but smaller, mostly negative trends in northeastern Canada (Fig. 2b).

River discharge
Although model-simulated discharges have previously been compared to measurements from major Arctic river watersheds under different climate conditions (Park et al. 2016(Park et al. , 2017;;Hiyama et al. 2023), the simulated discharges for the six major Arctic river watersheds in this study were compared against the measurements for model validation (Additional file 1: Fig. S1).The model generally well simulated the seasonal variation in discharge observed at the mouth of each watershed, particularly for the hydrograph of the spring peak discharge, fractioning the larger amount of the annual discharge.However, the model considerably overestimated the volume of the spring peak discharge in the Ob and Yenisey rivers.
The water flow in the upstream regions of the two river watersheds was disturbed by anthropogenic regulation from several artificial reservoirs (McClelland et al. 2004), for which a deficiency exists in the model representation, as mentioned in Sect.2.1, which likely resulted in the overestimations.In contrast, the model underestimated the observed summer discharge of the Yukon River, for which the upstream region is covered by a mountain glacier and where melted glacier water accounts for a considerable fraction of summer discharge (Chesnokova et al. 2020).This glacial process was not considered in the model simulation; thus, the deficiency was attributed to the bias in the simulated discharge.The degree of overestimation or underestimation of the discharge of the three rivers (i.e., Ob, Yenisey, and Yukon) was determined by the lowest Nash-Sutcliffe coefficients, which were compared to the coefficients with values of > 0.65 in the simulations of the remaining watersheds (i.e., Lena, Kolyma, and Mackenzie).
As indicated in Sect.2.1, surface and subsurface runoffs simulated by CHANGE were stored to individual reservoirs of the TRIP2 river channel then flowed down along a prescribed channel network.Figure 3 illustrates the seasonal variations, averaged over the study period, of the simulated daily runoff (Fig. 3a) flowing from the circumpolar Arctic river basins (Fig. 1) into the Arctic Ocean, and the corresponding seasonal trends (Fig. 3b).The total discharge (sum of the surface and subsurface runoffs) followed a typical hydrograph, with a marked seasonal peak in spring and early summer, followed by a quick decrease, then a smaller secondary peak in autumn (Fig. 3a).Peak discharge originated nearly exclusively from surface runoff, which accounted for 90% of the April-July total discharge.Conversely, subsurface runoff exhibited a different seasonality, increasing from August to its peak value in October, and represented the primary contribution to river discharge in the Arctic cold season (September-May).
Daily surface runoff, the primary contributor to total summer discharge in the Arctic, exhibited sharply positive than negative trends in the early summer season (May-June).Surface runoff greatly increased in spring, with a maximum value of 0.5 km 3 year −1 in late May (Fig. 3b), likely because of the T a increase.Increasing T a resulted in earlier spring snowmelt, with faster inflow of snowmelt into the circumpolar Arctic river discharge (Hiyama et al. 2023).This modification of the snow cycle seasonality has already been identified in observational data from the Eurasian and North American river basins (Mudryk et al. 2017(Mudryk et al. , 2021)).Then, the simulated surface runoff trend became negative, with a minimum value of − 0.2 km 3 year −1 at the end of June (Fig. 3b).We attribute this negative trend to a decrease in soil moisture caused by enhanced evapotranspiration, directly because of the increasing T a values.Evapotranspiration presented a trend of 0.8 km 3 year −1 larger than that of precipitation in June over the circumpolar Arctic river basin (Additional file 1: Fig. S2).Like its seasonal variations, the total river discharge trends were also determined primarily by surface runoff in spring and summer (Fig. 3b).In contrast, subsurface runoff exhibited a positive trend during the Arctic cold season, particularly after October (Fig. 3b), with a maximum value of 0.05 km 3 year −1 .The subsurface runoff increase likely resulted from large positive P r trends in summer and autumn (Additional file 1: Fig. S3).In previous studies, we had also discussed a groundwater (subsurface) contribution increase to the anomalous river discharge during the Arctic cold season, induced by the deepening of the active layer under warming climate conditions (Hiyama et al. 2023).
Annual mean discharge from the circumpolar Arctic river basins into the Arctic Ocean was 5 273 km 3 , increasing continuously at an annual mean rate of 12.2 km 3 year −1 (Fig. 4) during the study period.Corresponding increasing trends were also derived for the surface and subsurface runoffs.Interestingly, although surface runoff represented 71% of the total river discharge (notably more than subsurface runoff ), we calculated a smaller positive trend for surface runoff than for subsurface runoff, primarily explained by the markedly different rates of increase for P s and P r (Fig. 2c, see also next section).
Additional analysis for river discharge was conducted, breaking down the pan-Arctic scale to the watershed scale.Targeted to the major six Arctic river watersheds, the separated annual discharges were averaged for the study period, and their trends were analyzed (Table 2).As identified at the Arctic scale, the largest contribution of surface runoff to total river discharge, fractioned from 62 to 79% over the six river basins, was also consistent at the watershed scale.However, different to the discharge of the other five rivers that exhibited a trend of increase, the Yukon River discharge revealed a negative trend (− 0.22 km 3 year −1 ), although it was statistically insignificant.The negative trend in Yukon River discharge was likely influenced by the reduction in surface runoff discharge (− 0.36 km 3 year −1 ) linked to the decrease in snow water (Fig. 2a), as observed by satellite Fig. 3 Seasonal variability of the discharge volume and trends.Daily mean values calculated for the reference simulation over the study period are indicated in the study area for a the total, surface, and subsurface runoff discharges from the circumpolar Arctic river basins, and b the corresponding annual trends.Discharges associated with surface and subsurface runoffs simulated by CHANGE were calculated separately by TRIP2 along the prescribed river route network, then combined to derive the total discharge Fig. 4 Interannual discharge variability.Annual mean discharge values were calculated in the study area for the surface and subsurface runoffs routed separately in TRIP2, then combined to derive temporal variations of the total discharge.Corresponding trends, calculated over the study period, are also indicated remote sensing (Pulliainen et al. 2020).Additionally, the absolute value of the trend in subsurface runoff of the Lena River was larger than that of the surface runoff, but it deviated from the results derived for the other rivers.As previously documented by Hiyama et al. (2023), the increase in subsurface runoff of the Lena River is attributable to the influence of increased summer rainfall and permafrost thawing.Comparison of the discharge budget between the six river basins reveals regionally different changes in seasonal processes under the background of climate change.

Impact of snowmelt on the circumpolar Arctic river discharge
To evaluate the sensitivity of the surface and subsurface runoffs to altered climate conditions during the Arctic cold season, particularly to snowmelt variations, we calculated discharge anomalies by subtracting the reference simulation from the treated simulations (with T a and P as treatment variables, Table 1).Figure 5 illustrates the seasonal anomalous discharge variations for each simulation.In the three P-treated simulations (i.e., EX PR_SON , EX PR_DJF , and EX PR_MAM ), larger P values during the cold season, inducing larger P s , clearly resulted in larger surface runoff in spring and early summer, because of the increased snowmelt relative to the reference simulation (Fig. 5a).The largest anomalous surface runoff among the three P-treated simulations resulted from the EX PR_DJF simulation, while that of the EX PR_SON simulation was lowest.The P increase in the mid-winter season with the coldest temperature was linked to larger snowpack and then higher snowmelt in spring.The variance of the seasonal P increments resulted in the discrepancy in annual river discharge between the three P-treated experiments, where the anomaly of the amount of evapotranspiration is negligible.Conversely, higher T a values during the cold season mainly induced a decrease in the April-July surface runoff, as shown by the EX TA_SON and EX TA_DJF simulations.Under warmer conditions, the contribution of P s to P was reduced, resulting in a thinner snow cover during the Arctic cold season and therefore less snowmelt in spring.The negative effect of T a on snowmelt was further evident in the EX TA_SON simulation with a low P s rate, relative to that in the EX TA_DJF simulation.However, the seasonality of the anomalous surface runoff discharge in the EX TA_MAM simulation (treated T a values in spring) was inconsistent with the other T a -treated simulations (EX TA_SON and EX TA_DJF ); the temperature increase first caused an earlier (April-May) snowmelt-induced positive surface runoff anomaly, and then a strong negative anomaly (mid-June) that was probably the result of soil Fig. 5 Seasonal variability of anomalous discharges in the sensitivity simulations.The anomalies represent discharge differences associated with the a surface and b subsurface runoffs between the treated simulations (with increases in surface air temperature, precipitation, or both simultaneously) and the reference simulation (Fig. 3) Park et al. Progress in Earth and Planetary Science (2024) 11:13 dryness caused by temperature-enhanced evaporation from the soil surface.A similar pattern was identified in the EX TA_PR_MAM simulation, in which the combined T a and P increases resulted in larger early-spring surface runoff than in the EX TA_MAM simulation, mitigating the negative influence of T a on the surface runoff.The positive influence of P on the surface runoff in early summer was also identified in the other T a -and P-treated simulations (EX TA_PR_SON and EX TA_PR_DJF ), in which the anomalous surface runoff exhibited similar seasonality to that in the three P-treated simulations (i.e., EX PR_SON , EX PR_DJF , and EX PR_MAM ).
Our results also indicated significant sensitivity of the subsurface runoff to cold season climate variations when P increased in autumn (i.e., simulations EX PR_ SON and EX TA_PR_SON , Fig. 5b).Evapotranspiration in autumn, when atmospheric demand is lowest, is typically extremely low (Park et al. 2008;Kim et al. 2023).Under such conditions, excess precipitation effectively contributes to soil wetting, thus enhancing subsurface runoff.In contrast, the higher temperature in autumn increases evapotranspiration and results in negative anomalous subsurface runoff (i.e., the EX TA_SON simulation).In the remaining seven simulations other than the EX PR_SON and EX TA_PR_SON simulations, subsurface runoff sensitivity to cold season climate variations was moderate, with anomalous discharge values consistently within ± 20 km 3 (Fig. 5b).The thicker snow cover caused by a winter increase in P produces higher insulation, thus a deeper active layer (Park et al. 2015), in which more thawed soil water is available for subsurface runoff.However, the dependence of summer subsurface runoff on the winter snow conditions was limited, as identified by the EX PR_DJF and EX PR_MAM simulations (Fig. 5b).This suggests that subsurface runoff depended primarily on summer and autumn rainfall, as reported in a previous study (Hiyama et al. 2023).Meanwhile, the T a increase in the cold season resulted in slightly negative anomalous subsurface runoff, particularly in late summer and autumn (i.e., the EX TA_ SON and EX TA_MAM simulations), which was probably due to the influence of higher evapotranspiration-induced soil dryness.
To investigate our hypothesis on the possible relationship between the surface/subsurface runoffs and P s /P r , we compared the anomalous annual mean surface and subsurface runoffs with P s and P r anomalies, respectively, for all simulations (Fig. 6).The impact of snowmelt on the runoff anomalies persisted until summer and was not limited to spring (Fig. 5); therefore, the anomalous annual mean runoff values were used for comparison purposes.The comparison indicated a strong, statistically significant linear correlation (correlation coefficient r = 0.98, p < 0.001) between P s and surface runoff.On the basis of this linear relationship, we estimated directly the contribution of P s to the surface runoff: a P s increase of 10 mm, equivalent to a water volume of 211 km 3 , would account for 73.6% of the surface runoff discharge anomaly (286.8 km 3 ; Fig. 6).Moreover, a T a increase of 2 °C reduced P s by (3-13) mm, corresponding to surface runoff decreases of (89-240) km 3 .The second comparison between subsurface runoff and P r anomalies also yielded a statistically significant linear correlation, but with more data scatter and a lower correlation coefficient (r = 0.75, p < 0.01).Thus, the contribution of P r to subsurface runoff was more difficult to evaluate.However, from simulations EX PR_SON and EX TA_PR_SON , we estimated that P r increases of (20-38) mm in autumn could induce a subsurface runoff increase of approximately 328 km 3 .Finally, the total anomalous circumpolar Arctic river discharge was significantly correlated to the total annual P anomaly (Fig. 6e), suggesting that P was the main contributor to the yearly total circumpolar Arctic river discharge.Additional comparisons showed that there were no significant correlations between P s and subsurface runoff (Fig. 6c) or P r and surface runoff (Fig. 6b).
We also investigated the spatial distribution of the surface and subsurface runoffs to determine the main areas of circumpolar Arctic river discharge.For this purpose, we calculated correlations between annual mean P s and P r and CHANGE-simulated surface and subsurface runoffs, respectively, in each grid cell for all simulations (Fig. 7).The surface runoff distribution yielded statistically significant correlations with P s (r > 0.6, p < 0.05) in all Arctic river basin grid cells (Fig. 7a), implying that snowmelt was the largest contributor to surface runoff.As snow melts in spring, soil surface layers remain initially frozen (Park et al. 2015), effectively generating mainly surface runoff from the snow meltwater by preventing soil infiltration.Significant correlations (p < 0.05) were also calculated between subsurface runoff and P r , but they were regionally localized, mainly at the southern permafrost boundary and in non-permafrost circumpolar regions of the study area (Fig. 7b), in good correspondence with the increasing P r trend regions of Fig. 2b.Interestingly, these regions were also characterized by large standard deviations (over all simulations) associated with active layer thickness and soil moisture in the upper layers (Fig. 8), indicating larger variability of the diagnostic variables in the treated simulations.When the active layer deepens under the influence of stronger snow insulation or higher temperatures, the soil water storage capacity also increases (Suzuki et al. 2021).

Fig.
Relationship between anomalous discharge components and anomalous precipitation variables.Integrated annual anomalous values, calculated in each treated simulation, are plotted for the discharge components (surface and subsurface runoff, total discharge) against the precipitation variables (snowfall, rainfall, and total precipitation) to illustrate the relationships between a surface runoff-snowfall; b surface runoff-rainfall; c subsurface runoff-snowfall; d subsurface runoff-rainfall; and e total discharge-total precipitation.Anomalies represent annual mean differences between each treated simulation and the reference simulation (Fig. 5).Values for the precipitation variables are calculated as area-weighted averages in each grid cell of the circumpolar Arctic river basins.Numbers in symbols represent experimental scenarios: 1 EX TA_SON , 2 EX PR_SON , 3 EX TA_PR_SON , 4 EX TA_DJF , 5 EX PR_DJF , 6 EX TA_PR_DJF , 7 EX TA_MAM , 8 EX PR_MAM , and 9 EX TA_PR_MAM This study investigated the influence of snowmelt on circumpolar Arctic river discharge, on the basis of treated model simulations that incorporated climate modifications, with T a and P as treatment variables, in the Arctic cold season.Snowmelt in spring generally separated into two flows: surface runoff and soil infiltration.The onset of snowmelt occurred earlier than soil thawing (Park et al. 2015).Frozen soil limited snow meltwater infiltration, therefore meltwater mainly generated surface runoff, particularly at the initial stage of thawing.Surface runoff accounted for 90% or more of the peak Arctic river discharge in spring and early summer (Fig. 3a).The strong relationship between snowmelt and surface runoff was demonstrated by high correlations in all grid cells of the Arctic river basins (Fig. 7).This suggested that P s was the primary contributor to surface runoff change.We also calculated a small, but not statistically significant, increasing P s trend over the study period (0.03 mm year −1 , p < 0.41; Fig. 2c), insufficient to explain the strong surface runoff increasing trend (Fig. 4).The missing contribution was likely provided by P r increases in autumn of the previous year.Part of the autumn rainfall, captured and frozen in winter in upper soil layers, thawed during the following spring and mixed with infiltrated snow meltwater, thereby also contributing to surface runoff.This is commonly referred to as the "memory effect" of hydrological processes in the Arctic circumpolar regions (Park et al. 2021;Hiyama et al. 2023).
Precipitation over the circumpolar Arctic river basins exhibited an increasing trend over the study period (Fig. 2c) under warming climate conditions.Variations in T a critically affect the partitioning of P into P s and P r , with a stronger influence in spring and autumn when transitions between freezing and thawing occur.However, the measured P s and P r trends were negative and positive, respectively, in spring and autumn (Additional file 1: Fig. S2).The P r increase in spring possibly mitigated the surface runoff decreasing trend induced by a P s decrease.Conversely, decreasing P s values in autumn, when snow starts accumulating, corresponded to lower winter SWE values, thus to a lower contribution to spring surface runoff.Previous climate model studies have also projected a P s decrease in the Arctic circumpolar regions under warming climate conditions (Bintanja and Selten 2014), with reduced snow cover and shorter snowing period (Mudryk et al. 2021;Robinson 2021).This implies that the contribution of snowmelt to runoff generation will decrease under warming climate conditions, as also demonstrated by our T a -treated simulations (Fig. 6a).In the present-day climate simulation, summer P s exhibited a decreasing trend (Additional file 1: Fig. S2), likely inducing a reduction in the mountainous snowpack.However, it is difficult to directly attribute the snow loss in mountainous regions to a decreasing river flow, because of the replenishment effect induced on P r by a T a increase.Assessing these warming-induced effects is also complicated by competing factors such as higher evapotranspiration, altered vegetation composition, and changes in wildfire propagation.However, a warmer climate will demonstrably alter the dependence of runoff on snowmelt in the Arctic.Similar events have already occurred on Earth, as indicated by historical trends established at the global scale in mountainous regions (Huss et al. 2017;Immerzeel et al. 2020;Siirila-Woodburn et al. 2021).
To analyze the sensitivity of the circumpolar Arctic river discharge to snowmelt, we conducted model simulations in the Arctic cold season, with T a and P as treatment variables.The treated variables were sequentially updated following modifications of other meteorological variables.For example, higher T a altered specific humidity, thereby inducing a P increase that, in turn, resulted in a lower T a value.The simple treatment implemented in our simulations only altered the diagnostic variables but excluded internal feedbacks between meteorological variables.This simplification causes uncertainties on the simulated results: an increase in T a yields a higher saturated vapor pressure, thus a larger vapor pressure deficit, leading to an overestimation of snow surface sublimation, thereby introducing errors in SWE calculations.Furthermore, in real environmental conditions, frozen and thawed soil often coexist at the regional scale, a likely characteristic of transitional landscapes at the boundaries between continuous and discontinuous (sporadic or isolated) permafrost regions (Kim et al. 2012).In our model, the spatial resolution is coarse (0.5° × 0.5°) and simulated soil in each horizontal grid cell is assumed completely frozen or fully thawed, with direct implications for surface runoff and soil infiltration.Therefore, the limited model configuration cannot fully reproduce hydrological processes at the soil surface in transitional boundary regions, causing errors on the simulated surface runoff.For example, our simulation results indicated that the P s increase accounted for 73.6% of the total Arctic river discharge variations (Fig. 6), larger than the 60% contribution previously measured in the Kolyma and Yukon river basins (Welp et al. 2005;Lachniet et al. 2016).These Fig. 8 Spatial distribution of standard deviations for a active layer thickness and b surface soil moisture.Surface soil moisture is defined as the water content in a soil column extending from the surface to a depth of 1.5 m.The standard deviations, calculated from the anomalous values representing the difference between the treated simulations and the reference simulation, illustrate the sensitivity of the active layer and the soil moisture to climatic variations during the Arctic cold season Park et al. Progress in Earth and Planetary Science (2024) 11:13 discrepancies likely result from limitations in the experimental design and model representation.Snowmelt in spring linked to snow meltwater, and soil surface water attributable to P r in autumn of the previous year, are contributors to surface runoff.The estimated value of 73.6% could include a partial contribution of P r , indicating a limitation of the current model in that it cannot separate the individual contributions of snow meltwater and rain water to surface runoff.However, results obtained in this study were consistent with observed historical snow loss trends induced by climate warming, despite the simple simulation settings.Therefore, we believe that our simulation results are useful to quantify regional water budgets and analyze the influence of snow conditions on water discharges in circumpolar Arctic river basins under warming climate conditions.

Conclusions
The total discharge inflow from circumpolar Arctic river basins into the Arctic Ocean has increased continuously in the last decades.This study quantitatively evaluated the contribution of snowmelt, sensitive to climate variations, to the observed increase in Arctic river discharge.Using model simulation results for 1979-2018, we established that snowmelt, through surface runoff, was a major contributor to both the spring peak and total annual river discharge.Correlation between snowmelt and surface runoff was strong in nearly all grid cells of the circumpolar Arctic river basins.Treated simulations, initialized with fixed increases in T a , P, or both to simulate climate warming variations, indicated that an increase in T a enhanced snowmelt, thereby reducing surface runoff and river discharge, but also induced larger rainfall and subsurface runoff, particularly in the autumn season.These results synthetically demonstrate that climate warming will limit the contribution of snowmelt to the total Arctic river discharge while increasing the influence of rainfall on the discharge.However, the simplified model assumption did not account for feedbacks between meteorological variables; therefore, our results also include the uncertainty generated by this simplification.Using climate model projections to represent the near-future warming climate conditions might be beneficial to reduce the uncertainty caused by the model settings' simplification.One of the biggest problems in the field of cold land hydrology is to determine the changes in seasonal processes under the background of a warming climate.
Therefore, improving the model to separate the variations in source waters attributable to the seasonal hydrology is a priority for our future work.

Fig. 2
Fig. 2 Spatial distribution and interannual variability of precipitation (snowfall and rainfall).Spatial distributions of the mean annual trends during the study period are represented in the study area (illustrated in Fig. 1) for a snowfall and b rainfall, in which hatching indicates statistically significant trends at the 90% confidence level.Panel c corresponding interannual variability for total precipitation, snowfall, and rainfall.The mean annual trends (values indicated in the figure with the associated p values) were calculated as area-weighted averages in each grid cell of the Arctic circumpolar river basins

Fig. 7
Fig. 7 Spatial distribution of correlations between runoff and precipitation.Correlation coefficient values are represented for a surface runoffsnowfall and b subsurface runoff-rainfall correlations.The correlation coefficient was calculated in each grid cell of the study area from the annual mean values in each treated simulation (nine values per grid cell).It represents the influence of snowfall on surface runoff (Panel a) and of rainfall on subsurface runoff (Panel b).Correlations are only plotted in the grid cells where statistical significance exceeds 5%

Table 1
Model settings for the sensitivity simulations *The differences adopted in this study for the climatic treatment variables (second and third columns) represent increases relative to the reference simulation **n/a not

Table 2
Summary of discharge budget for the six major Arctic river basins