Model Diagnostic Analysis in a Cold Basin Influenced by Frozen Soils With the Aid of Stable Isotope

Understanding the hydrological processes on the Tibetan Plateau (TP) under climate change is an important scientific question. The frequent multiphase transfer exacerbates the complexity of hydrological processes on the TP, which brings equifinality problem to hydrological models and causes large uncertainties in quantifying the contributions of runoff components. Tracer‐aided hydrological models are helpful for improving model performances and have been adopted in cryospheric regions, but the influence of frozen soil has yet to be considered. This study adopted the Tracer‐aided Tsinghua Representative Elementary Watershed model (THREW‐T) in a typical cold basin with widespread frozen soil on the TP. The model structure was diagnosed with isotope by identifying the influences of frozen soil. A simplified catchment‐scale frozen soil module was incorporated into the model. Results showed that: (a) The THREW‐T model cannot simultaneously simulate baseflow and stream water isotope well. The imbalance of simulations on two objectives could be attributed to the influence of frozen soil, resulting in seasonal variation of soil‐related parameters, which was not considered in the model. (b) Incorporating the frozen soil module significantly improved the balance of baseflow and isotope simulation, simultaneously producing low baseflow and high contribution of subsurface runoff during wet seasons. (c) The frozen soil had little influence on the annual streamflow, but changed the runoff seasonality by reducing baseflow during dry seasons and increasing subsurface runoff during wet seasons. The frozen soil module was still simplified, and further work is needed to improve the physical representation of soil freeze‐thaw process. This study highlights the value of tracer‐aided hydrological modeling method on diagnosing model structure by identifying the influence of specific processes such as frozen soil.


Introduction
The Tibetan Plateau (TP), termed as the "Asian Water Tower" and the "Third Pole" (Su et al., 2023), is the source of several Asian major rivers such as Indus, Brahmaputra, Mekong and Yangtze, playing a crucial role in providing essential freshwater resources for downstream livelihoods, ecosystems and agricultural irrigation (Immerzeel et al., 2010;Zhang et al., 2013).Meanwhile, the air temperature has warmed much faster than the global average, making hydrological processes on the TP response sensitively to climate change (Li et al., 2021).
A comprehensive understanding of the hydrological processes on the TP holds significant importance not only for effective catchment management but also for its broader scientific relevance in the context of climate change.
As a typical mountainous cryosphere, the TP exhibits widespread distribution of frozen water, leading to frequent water phase transfer (Li et al., 2019;Yao et al., 2022).This results in the complexity of hydrological processes on the TP, characterized by the diverse runoff components (Cui et al., 2023).The extension-induced rift systems on the TP may transport substantial quantities of groundwater to rivers and induce inter-basin groundwater flow (Tan et al., 2021;Yong et al., 2021), further exacerbating the complexity of hydrological processes.Accurately estimating the contributions of runoff components to streamflow is crucial for a reliable projection of future runoff change (Cui et al., 2023;Lutz et al., 2014;Xu et al., 2019).However, this is a challenging task (Weiler et al., 2018).The compensatory effects of the runoff induced by different water sources are often inadequately constrained in the hydrological models (Duethmann et al., 2015), resulting in equifinality problem (Gupta et al., 2008) and large variations in the runoff component estimations.Taking the source region of Brahmaputra River as an example, the contribution of glacier meltwater to annual runoff reported in existing publications varies widely, ranging from 3.5% to 29% (Boral & Sen, 2020;Wang et al., 2021).Moreover, the reported contribution of groundwater runoff even has a larger uncertainty, ranging from 4% to 80% (Liu et al., 2020;Wang et al., 2021).Consequently, reducing the equifinality of hydrological models and accurately estimating the contributions of runoff components stand as significant scientific challenges in understanding the hydrological processes on the TP.Apart from the above processes, the soil freeze-thaw process also influenced the hydrological processes on the TP significantly.Climate warming has caused widespread frozen soil thawing (Schuur & Mack, 2018), which have severe consequences for the ecological and hydrological systems (O'Neill et al., 2020;Wang & Gao, 2022).The presence of both liquid water and ice in frozen soil significantly modifies the hydraulic and thermal properties of the soil (Decharme et al., 2016), which has significant influences on the water and energy exchange in the atmosphere-land system (Stevens et al., 2007).The freeze-thaw cycle of frozen soil modifies the soil permeability and changes the magnitudes and temporal patterns of runoff (Qin et al., 2016;Wang et al., 2017).The frozen soil processes consequently increase the hydrological complexity in cold regions, and form an important uncertainty source of hydrological model and runoff component estimation.Frozen soil modules have been developed for hydrological and earth system models to simulate the coupled heat and water transfer processes, which could be grouped into two categories: analytical and numerical approaches (Gao et al., 2022).The Stefan equation is a widely used analytical method to calculate the frozen soil depth and the impacts of frozen soil on hydrology based on ground temperature and thermodynamic properties of soil (Fabre et al., 2017;Votyakova, 1976;Xie & Gough, 2013).Such kind of method allows for a simplified description on the frozen soil, with the advantage of requiring less data and the simple calculation process, but it cannot produce the detailed soil characteristics over the soil profile.A typical example of frozen soil hydrological model is the FLEX-Topo-FS model developed by Gao et al. (2022), which calculated the frozen soil depth and the volume of frozen groundwater using the Stefan equation, but the impact of frozen soil on soil hydraulic properties were not considered.Numerical solution schemes such as finite difference are typically adopted in earth system and land surface models to produce detailed simulations on the water content and temperature over the soil profile and the changes in frost-thaw fronts (e.g., Gao et al., 2019;Li et al., 2023;Liu et al., 2013;Zheng et al., 2019).While numerical methods offer detailed process descriptions, they require extensive data and high-resolution spatio-temporal discretization, making them unsuitable for conceptual or semi-distributed hydrological models.
Tracer-aided hydrological modeling method that integrates environmental tracers (such as stable oxygen isotope 18 O) have been proven useful for improving hydrological simulations (Birkel & Soulsby, 2015).This is generally achieved by either constraining the parameter uncertainties (He et al., 2019;McGuire et al., 2007) or diagnosing the model structures (Birkel et al., 2011;Son & Sivapalan, 2007).The tracer-aided models characterize the water storage, mixture and transportation processes, as reflected by the "velocity," by simulating the tracer movement related to water fluxes (McDonnell & Beven, 2014).Different runoff components generally have different isotope signatures because of the varying mixture and phase change processes.For instance, the isotope composition of subsurface runoff with longer travel time tends to be more damped compared to the quick surface runoff (McGuire & McDonnell, 2006), while the isotope composition of snow/glacier meltwater is usually more depleted compared to rainfall water due to the isotopic fractionation (Xi, 2014).The stream water, as the output of catchment hydrological processes, is the mixture of multiple components, so the difference in isotope composition among various water bodies can provide additional information to determine the contribution proportion of each component.Tracer-aided models consequently show significant value in regulating the internal hydrological processes and constraining the water apportionments among various components (Delavau et al., 2017;He et al., 2019;Nan, Tian, et al., 2021).With the deeper understanding of the influence of cryospheric processes on isotope (Nan et al., 2022;Pu et al., 2020) and the development of isotope-enabled climate circulation models (Noone & Sturm, 2010), the tracer-aided models have been increasingly applied in cold mountains basins with different drainage areas, where the hydrological processes are typically more complex (Ala-aho et al., 2017;He et al., 2019;Nan, He, et al., 2021, Nan, Tian, et al., 2021).
Although the tracer-aided modeling method beard the potential to diagnose model structures, most of recent tracer-aided hydrological modeling studies focused on improving model performance by constraining parameter uncertainties through the multi-objective calibration strategy (Holmes et al., 2023;Stadnyk & Holmes, 2023) or attempted to gain deeper insight into the runoff generation mechanisms (Wang et al., 2023).Very few studies aimed to diagnose the model structure and improve the physical representation.For instance, as one of the few model diagnosis studies, Birkel et al. (2011) utilized isotopes to identify the effects of passive storage and groundwater return flow processes, and incorporated the modules related to these processes to improve model performance.The value of tracer-aided modeling method on improving model structure is poorly explored.Besides, although tracer-aided models have been applied successfully in cold mountainous catchments, the influence of frozen soil is yet to be addressed in the existing tracer-aided modeling studies, and its impact on the isotope composition of water is poorly understood.The soil freeze-thaw processes regulate the partitioning among different runoff generation pathways by influencing the soil properties such as hydraulic conductivity and water permeability (Gao et al., 2022;Song et al., 2022).Such impact is expected to be effectively identified by the isotope signature due to the significant difference in the isotopic composition of different runoff components, which is however poorly explored by now.In our previous study in the source region of Yangtze River (SRYR), a typical large basin on the TP with widespread frozen soil, we found that the streamflow during dry seasons was extremely low, difficult to simulate when other objectives including isotope were reproduced well, which might be the result of frozen soil (Nan et al., 2023).Motivated by the above background, we conducted a model diagnostic study in the source region of the Yangtze River adopting the Tracer-aided Tsinghua Representative Elementary Watershed model (THREW-T) developed by Nan, He, et al. (2021).The objective of this study is to test the value of tracer-aided modeling method on identifying the hydrological impact of frozen soil.This is achieved by testing the model performance on streamflow and isotope simulation, analyzing the reason for the imbalanced performance on multiple objectives and developing a frozen soil module incorporated into the THREW-T model to characterize the influence of soil freeze-thaw on hydrological processes.
The paper is organized as follows: Section 2 briefly introduces the study area and data.Section 3 describes the THREW-T model and the calibration variants.Section 4 presents the performance of original THREW-T model, and analyzes the reason for its bad performance.Section 5 proposes a frozen soil module and presents the performance of model incorporated with the module.Section 6 discusses the limitations and novelties of this study.The conclusions are finally drawn in Section 7.

Study Area and Data
The SRYR was selected as the study region of this paper.It is located in the hinterland of the TP, spanning from 32°N to 36°N and 90°E to 97°E (Figure 1).The SRYR basin has an elevation range of 3,600-6,300 m, with an average elevation of 4,750 m.The outlet station of SRYA basin is Zhimenda station, with a drainage area of approximately 138,000 km 2 .The average annual precipitation is approximately 440 mm, most occurring in summer months.The mean annual temperature is 5.1°C, which is the lowest among the major river basins on the TP (Cui et al., 2023), resulting in a widespread distribution of permafrost and snow.The permafrost covers approximately 70% of the total basin area (Shi et al., 2020;Zou et al., 2017), significantly influencing the hydrological processes in the SRYR basin (Gao et al., 2012;Li et al., 2020).The average snow cover area (SCA) is around 15% of the total basin, while glaciers only cover approximately 0.7% of the basin (Yao et al., 2014).
We compiled the data sets representing the topographic, meteorological and underlying conditions in the SRYR basin to establish the model.The 30 m digital elevation model (DEM) data were downloaded from the Geospatial Data Cloud website (https://www.gscloud.cn).The daily precipitation and temperature data were extracted from the 0.1°China Meteorological Forcing Dataset (CMFD, Yang & He, 2019).The daily potential evaporation data was extracted from the ECMWF Reanalysis v5 -Land (ERA5-Land, Munoz-Sabater et al., 2021) with a resolution of 9 km.The daily TP Snow Cover Extent (TPSCE) product (Chen et al., 2018) was adopted for the evaluation of snow simulation.The second glacier inventory data set of China (Liu, 2012) was adopted to identify the glacier cover regions.Vegetation situations were reflected through 8-day Leaf Area Index and monthly Normalized Difference Vegetation Index from MODIS satellite products MOD15A2H (Myneni et al., 2015) and MOD13A2 (Didan, 2015).Soil parameters were estimated from the Global high-resolution data set of soil hydraulic and thermal parameters (Dai et al., 2019).Daily streamflow data from 2010 to 2018 at Zhimenda station were collected to evaluate the streamflow simulation.
Stream water samples were collected weekly during the ablation period of 2016 and 2017 at Zhimenda and Tuotuohe stations (Figure 1b).These samples were analyzed for isotope compositions, and the basic characteristics of the samples, along with the isotope data of stream water, are detailed in Table 1.The simulated precipitation isotope by Scripps global spectral model incorporated with water isotopes (isoGSM, Yoshimura et al., 2008) were adopted as the isotopic input data for the tracer-aided model, to denote the spatiotemporal variation of the isotope composition in precipitation.The isoGSM product has been successfully utilized as precipitation isotope input data for isotopic-hydrological modeling across various large basins (Arciniega-Esparza et al., 2023;Nan et al., 2022Nan et al., , 2023)).Previous evaluations of the isoGSM based on the monthly precipitation isotope data at 19 stations over the TP (Figure 1a; Gao, 2020) indicated that while the isoGSM product effectively captured the seasonal fluctuation of precipitation δ 18 O, it exhibited limitations in accurately producing the isotopic variations of specific events or periods (Nan, He, et al., 2021, 2023).The isoGSM tended to overestimate the precipitation δ 18 O in low-latitude regions of TP and underestimate it in high-latitude regions.The regression equation between the isoGSM bias and geographical factors (including latitude, longitude and elevation) obtained by Nan et al. (2023) was adopted to correct the isoGSM.The spatial and temporal patterns of the precipitation δ 18 O derived from the corrected isoGSM are shown in Figure 2.More details of the bias characteristic and correction equation of isoGSM can be found in Nan et al. (2023).

Description of the THREW-T Model
In this study, we adopted a distributed tracer-aided hydrological model, Tsinghua Representative Elementary Watershed-Tracer-aided version (THREW-T) model, developed by Nan, Tian, et al. (2021) to simulate the hydrological processes and isotopic variations.The model uses the REW method (Reggiani et al., 1999) to divide the catchment into several REWs based on the DEM data, and each REW is further divided into surface layer and subsurface layer (Tian et al., 2006).The surface layer consists of six subzones with different underlying types: vegetation zone (v-zone), bare soil zone (b-zone), main channel reach zone (r-zone), sub-stream network zone (tzone), snow zone (n-zone) and glacier zone (g-zone).The subsurface layer includes two subzones according to soil saturation conditions and soil depths, namely the unsaturated zone (u-zone) and the saturated zone (s-zone).
Figure 3 illustrates the schematic representation of the THREW-T model structure.In this study, the SRYR basin was divided into 178 REWs,similarly with Nan et al. (2023).More details of the THREW-T model can be found in Tian et al. (2006) and Nan, Tian, et al. (2021).
A cryospheric module was developed for application in cold regions (Nan, Tian, et al., 2021).Precipitation was partitioned into liquid (rainfall) and solid (snowfall) based on a temperature threshold set as 0°C.The degree day factor method was adopted to calculate the meltwater from n-zone and g-zone.The snow water equivalent of each REW was updated based on the snowpack mass balance, and the SCA was determined by the snow cover depletion curve.Considering the extremely small glacier cover area in the SRYR basin, the area of glacier covered region was assumed constant for simplification.The tracer module was integrated to simulate the isotope composition in multiple water bodies.The Rayleigh equation was used to characterize the isotope fractionation during evaporation and snow melting processes, similar with He et al. (2019).A complete-mixture assumption within each simulation units was adopted to calculate the mass and concentration of tracers, similar with most of distributed tracer-aided hydrological models (e.g., Ala-aho et al., 2017;Delavau et al., 2017;He et al., 2019).The isotope composition of glacier meltwater was assumed as a constant value, lower than the local average precipitation isotope composition by an offset parameter (Nan et al., 2022).The isotope fluxes among different simulation units and the isotope compositions were determined based on the mass balance equations by the water flux variables, which were already calculated in the original hydrological module of THREW-T model.
The THREW-T model defined the runoff components based on the runoff generation flow-pathway as reviewed by He et al. (2021).The total runoff was divided into surface and subsurface runoff (baseflow).The surface runoff included the rainfall, snowmelt and glacier melt occurring in the saturation excess runoff area, which was determined by Equation 1 (Tian et al., 2008).The subsurface runoff is the groundwater outflow, which was calculated by Equation 2 (Tian et al., 2012): (1) where, A s is the saturation excess runoff area.WM and B are the shape coefficient and tension water storage capacity, respectively.y s and y u are the water depth in s-zone and u-zone, respectively.s u is the saturation degree of u-zone.ε u and ε s are the volumetric soil moistures in u-and s-zone, respectively.R g is the groundwater outflow rate.S is the topographic slope.KKD and KKA are the linear and exponential coefficients, respectively.K s s is the saturated hydraulic conductivity.Z is the total soil depth.Among these parameters, WM, B, KKD, and KKA are determined by calibration.y u , y s and s u are state variables which are updated by the water mass balance equations.The remaining parameters are determined according to the soil and topography data sets.

Model Calibration
The physical meanings of the calibrated parameters are listed in Table 2. Two variants were designed for the calibration procedure: (a) double-objective calibration toward discharge and SCA (NoFS-DS calibration variant), and (b) triple-objective calibration toward discharge, SCA and stream water δ 18 O (NoFS-DSI calibration variant).Considering the data availability, data sets used for calibration included: the discharge data at Zhimenda station during 2010-2018, SCA extracted from the TPSCE data set for the SRYR basin during 2010-2016, and the stream water δ 18 O data at Zhimenda and Tuotuohe stations during the wet season in 2016 and 2017.The Nash-Sutcliffe efficiency coefficient (NSE) and the logarithmic Nash-Sutcliffe efficiency coefficient (lnNSE) were used for discharge calibration to evaluate the simulation of both high and low flow processes.For SCA and isotope calibration, the root mean square error (RMSE) and NSE were adopted, respectively.These evaluation metrics were calculated by Equations 3-6.In DS and DSI calibration variants, the objective functions of the calibration procedure were set as Equations 7 and 8, respectively.
(5) where, n, m and k are the numbers of discharge, SCA and isotope measurements used for calibration.Subscripts "o" and "s" refer to the observed and simulated variables, respectively.
The Python Surrogate Optimization Toolbox (pySOT, Eriksson et al., 2017) was adopted for calibration.It adopted the symmetric Latin hypercube design strategy to generate parameter values, allowing an arbitrary number of design points.Each pySOT running would stop after a maximum time of model running was reached, which was set as 3,000 in this study.To generate sufficient parameter samples for capturing the optimal fronts, we repeated the pySOT for 100 times in each calibration variant.A final parameter set with the best objective function (Equations 7 and 8) was selected from the calibrated parameter sets.

Model Performance in Different Calibration Variants
Figure 4 shows the model performance obtained by NoFS-DS and NoFS-DSI calibration variants.When not calibrating isotope objective, the hydrograph was reproduced well in terms of both high and low flow processes (Figure 4a), as indicated by the high values of NSE and lnNSE.The variation of discharge was mostly captured by the model: There were generally more than one discharge peaks in each year, occurring during July-September; The discharge decreased sharply in October and remained extremely low during December-April.However, the simulated discharge from December to March was still overestimated, which could be about twice the observation.The SCA was simulated well by both two calibration variants, with very low RMSE (∼0.13 and ∼0.08 for daily and monthly SCA, respectively).The SCA rapidly decreased in May, remained almost zero in summer, and then increased quickly in September (Figures 4b and 4e), which was captured well in both variants.The performance of isotope simulation was extremely poor when not calibrating isotope, showing a much stronger variability than observation (Figure 4c), with a simulated standard deviation (1.51‰) higher than obversion (1.08‰), and a simulated average δ 18 O ( 12.28‰) lower than observation ( 11.31‰).When isotope was integrated into the calibration objective, its simulation was improved significantly, with NSE of 0.49.The stream water isotope variation during the ablation period was captured, characterized by an increase in June and early July, followed by a decrease in July and August, and subsequently another increase after September (Figure 4f).The simulated standard deviation (0.84‰) was slightly lower than the observation, and the simulated average δ 18 O ( 11.42‰) was very close to the observation.However, the discharge simulation got significantly worse, especially for the baseflow process, with an extremely low lnNSE of 0.44 (Figure 4d).The simulated discharge could be larger than 3.5 times the observation during December-February.Based on the above results, we can conclude that the current THREW-T model cannot simultaneously produce good simulation on the baseflow during dry seasons and stream water isotope during wet seasons.
The calibrated parameter values obtained by two calibration variants are listed in Table 3.The NoFS-DS calibration produced a higher B than NoFS-DSI, indicating a stronger spatial heterogeneity in the water storage capacity, which would result in a larger saturation excess area.Two calibration variants estimated similar KKA, but the KKD estimated by NoFS-DSI variant was significantly higher than that of NoFS-DS, indicating a larger groundwater outflow rate, which was consistent with the difference in baseflow simulation by two variants (Figures 4a and 4d).The melting rate simulated by NoFS-DSI variant was slower than NoFS-DS due to the lower DDF S /DDF G and higher T 0 , consistent with the simulated higher SCA during ablation season by NoFS-DSI (Figure 4e).Despite the large difference in parameters related to melting processes, the snow simulation produced by two variants were both acceptable, showing as the low RMSE.

The Reason for Poor Model Performance
The reason for the poor model ability on simultaneously simulating baseflow and stream water isotope was analyzed.The main issue with the isotope simulation was identified as a larger isotopic variability compared to the observations (Figure 4c).The isotopic variability of stream water reflected the water travel time and young water fraction within a catchment, which is influenced by water mixture processes (Jasechko et al., 2016;Kirchner, 2016;McGuire & McDonnell, 2006).To better understand the influence of runoff components on isotope simulation, the isotope signatures of different runoff components produced by the NoFS-DSI calibration variant are shown in Figure 5.The isotopic variability of surface runoff was quite strong because of its short travel time, but it was slightly damped compared to the precipitation input signature.This damping effect is a result of the water mixture occurring in the t-zone storage and the mixing of rainfall and snowmelt.On the contrary, the isotope signature of subsurface runoff was extremely steady, almost remaining constant, which was consistent with observation or simulation results of other studies (e.g., Birkel et al., 2011;He et al., 2019).The stream water was a mixture of these two components, and its isotope signature was within the range of two signatures (Figure 5).Consequently, the two components need to mix in a proper proportion to match the simulated isotope variability with the observation.
The contributions of subsurface runoff in the annual runoff estimated by NoFS-DS and NoFS-DSI calibration variants were 32.6% and 61.7%, respectively.According to Li et al. (2020) analyzing the runoff component contribution in the SRYR based on isotope data of around 1,000 water samples, the contribution of subsurface runoff was 31.52%-85.24%,varying by region and time.The proportion estimated by NoFS-DS calibration variant was at the lower end of this range, indicating an underestimation of subsurface runoff, which was consistent with the simulated stronger isotopic variability in Figure 4c.On the contrary, the NoFS-DSI calibration variant satisfied the isotope simulation by producing a large amount of subsurface runoff, but resulted in a high baseflow in dry season.Considering the widely distributed frozen soil in the SRYR basin (Li et al., 2020;Shi et al., 2020), the unbalanced simulation on baseflow and isotope could be attributed to the impact of frozen soil on the hydrological processes, which was not considered in the THREW-T model.The soils were generally frozen during dry seasons and thawing during wet seasons, which influenced the partitioning between surface and subsurface runoff in two ways.On the one hand, the frozen soil reduced the hydraulic conductivity, resulting in less subsurface runoff because of decreased groundwater outflow rate (Finney et al., 2012;Gao et al., 2018;Qi, Zhang, & Wang, 2019).On the other hand, the frozen soil reduced the soil water storage capacity, resulting in more surface runoff due to a larger saturation excess area (Bibi et al., 2019;Li et al., 2020).Likewise, frost thawing would have the opposite influence on hydrological processes in the wet season.The original THREW-T model did not consider the seasonal variation of these soil properties, leading to either consistently low or high subsurface runoff produced by NoFS-DS and NoFS-DSI calibration variants.Consequently, to achieve a simultaneous good performance in baseflow and isotope simulation, it is essential to simulate both low subsurface runoff in dry seasons and high subsurface runoff in wet seasons.Incorporating a frozen soil module into the model would be a feasible way to achieve this balance.

Description of the Frozen Soil Module
To characterize the hydrological impact of frozen soil, a correction factor (CF FS ) was introduced to establish the frozen soil module.CF FS was defined as the reduction coefficient of the soil-related parameters when frozen soil existed.The parameters WM and KKD in the THREW-T model were selected as soil-related parameters, because they reflected the soil storage and hydraulic conductivity properties respectively.They were multiplied by CF FS to calculate the parameter values influenced by frozen soil (Equations 9 and 10).
where, WM 0 and KKD 0 are the parameter values when there is no frozen soil, and WM FS and KKD FS are the values influenced by frozen soil.
According to the definition of CF FS , it should be related to the frozen soil conditions, and the value should be smaller when more frozen soils exist.Considering the general correlation between the freezing/thawing front depth and the square root of accumulated temperature (Ding et al., 2000;Hayashi et al., 2007;Qin et al., 2016), which was widely used in frozen soil simulation methods (e.g., Stefan equation, Votyakova, 1976), CF FS was assumed to be correlated with the accumulated temperature.Specifically, CF FS was set as 1 when there was no frozen soil (i.e., the positive accumulated temperature was higher than a certain value), and was assumed to be a minimum value CF m FS when the soil was adequately frozen (i.e., the negative accumulated temperature was lower than a certain value).CF FS was calculated by Equations 11-13 within each continuous thawing/freezing period with the same sign of daily temperature (i.e., positive or negative, respectively).
where, AT is the accumulated temperature, which is calculated from the beginning of each thawing/freezing period.ΔCF + FS and ΔCF FS are the changes of CF FS within each thawing and freezing period, respectively.AT + M and AT M are the maximum positive and negative accumulated temperature corresponding to the maximum range of CF FS .CF 0 is the initial CF FS at the beginning of the current thawing/freezing period, which is also the final CF FS at the end of the last freezing/thawing period.Three additional calibrated parameters AT + M , AT M , and CF m FS were introduced for frozen soil simulation.The THREW-T model incorporated with frozen soil module was calibrated toward all three objectives including discharge, SCA and stream water isotope (namely, FS-DSI calibration variant).The calibration procedure of FS-DSI variant was same with that of NoFS-DS and NoFS-DSI variants, as described in Section 3.2.
To better reflect the thermal condition of soil, the ground temperature rather than air temperature was adopted to calculate the accumulated temperature for frozen soil simulation.The REW scale mean annual ground temperature (MAGT) was extracted from the permafrost data sets over the Northern Hemisphere (NIEER, Ran et al., 2022).The REW scale MAGT was strongly correlated with the REW scale mean annual air temperature (MAT) extracted from the CMFD data set (Figure 6).This regression equation was adopted to calculate the daily MAGT as the input data for the frozen soil module.

Model Performance With the Frozen Soil Module
Figure 7 shows the simulations produced by the THREW-T model incorporated with the frozen soil module.The frozen soil module brought an overall improvement to the simulation on multiple objectives.Specifically, the simulation on baseflow process was improved significantly with an extremely high lnNSE of 0.9, higher than that produced by NoFS-DS calibration variant (0.77).The hydrographs of recession processes during October and November and baseflow during December-March were captured by the model almost perfectly (Figure 7a).The NSE of discharge simulation (0.81) was slightly lower than that of NoFS-DS (0.82), but the high flow processes were still reproduced well.Meanwhile, the simulations on SCA and isotope were slightly improved compared to NoFS-DS and NoFS-DSI calibration variants, showing as the lower RMSE for monthly SCA and higher NSE for isotope (Figures 7b and 7c).
The parameters obtained during the calibration procedure of NoFS-DSI and FS-DSI variants were saved to look into the relation among different simulation objectives.The envelope line of the relationship between NSE I and lnNSE D is shown in Figure 8.For the original THREW-T model (blue line in Figure 8), NSE I and lnNSE D decreased sharply with the increasing of another objective.When isotope simulation was at a high level with NSE I of ∼0.5, the baseflow simulation was poor with lnNSE D lower than 0.5.Similarly, when baseflow was simulated well with lnNSE D of ∼0.8, NSE I was lower than 0.2.On the contrary, incorporating the frozen soil module significantly improved the balance of isotope and baseflow simulation (red line in Figure 8).Good simulation on isotope and baseflow could be achieved simultaneously with NSE I and lnNSE D of 0.5 and 0.9, respectively.The contributions of subsurface runoff estimated by the model incorporated with frozen soil module were analyzed.The contribution of subsurface runoff in annual runoff was 59.0%, much higher than the estimation by NoFS-DS calibration variant (32.6%) and slightly lower than that of NoFS-DSI variant (61.7%).Figure 9 shows the monthly discharge of subsurface runoff and its contribution in total runoff.The FS-DSI variant (red line) estimated the discharge of subsurface runoff similarly to the NoFS-DS variant during dry seasons (blue line) and similarly to the NoFS-DSI variant during wet seasons (green line) simultaneously (Figure 9a).Meanwhile, the contribution ratio of subsurface runoff to total runoff during wet seasons estimated by FS-DSI calibration variant was also close to that of NoFS-DSI variant (Figure 9b).Consequently, the model incorporated with frozen soil module could produce good simulation on baseflow and isotope simultaneously.

Influences of Frozen Soil on Hydrological Processes
Figure 10 shows the time series of basin scale average daily ground temperature, accumulated ground temperature and CF FS estimated by the FS-DSI calibration variant.The MAGT was lower than 0°C (Figure 10a), so the yearly accumulated temperature kept decreasing (Figure 10b).However, the CF FS did not show interannual variability, because the calibrated AT + M was smaller than AT M , making CF FS rising back to ∼1 during the thawing season every year despite the smaller positive accumulated temperature (Figure 10c).In response to the seasonal temperature changes, CF FS exhibited distinct patterns: an increase during March-June, a relatively stable value of around 1 during June-September, a decrease during September-December, and a consistent presence around the minimum value (∼0.25) during December-March (Figure 10c).
To analyze the hydrological effects of frozen soil, we compared the simulations on multiple objectives in two scenarios.The model was driven by the parameters obtained by FS-DSI calibration variant to represent the condition with the influence of frozen soil (namely, FS scenario).In the scenario without the influence of frozen soil (namely, NoFS scenario), the model was also driven by the parameters obtained by FS-DSI calibration variant, but CF FS was fixed as the average value.Figure 11 shows the comparison of simulated discharge, stream water isotope and runoff component produced by two scenarios.Results indicated that whether considering the frozen soil or not had little influence on the total runoff, with mean annual discharge of 558.12 and 558.64 m 3 /s for the FS and NoFS scenarios, respectively.The runoff seasonality differed in two scenarios, with a higher discharge during dry season and a lower discharge during wet season in the NoFS scenario than in the FS scenario (Figure 11a).This was mainly a result of the higher groundwater outflow rate during dry seasons, and the lower rate during wet seasons in NoFS scenario.The difference in the seasonality of groundwater outflow rate also influenced the contribution of runoff components.Although the contribution of subsurface runoff in annual total runoff was similar in two scenarios (59.0% and 59.1% in the FS and NoFS scenarios), the seasonality of runoff component was different.Compared to the FS scenario, the amount and contribution ratio of subsurface runoff estimated by the NoFS scenario was higher during dry seasons, and lower during wet seasons (Figures 11c and 11d).As a result, the variation of stream water isotope during wet seasons was stronger in NoFS scenario (Figure 11b).In summary, the comparison indicated that the existence of frozen soil had little influence on the annual runoff, but significantly influenced the runoff seasonality, showing as a lower discharge in dry seasons and a higher discharge in wet seasons.The hydrological effect of frozen soil revealed by the THREW-T model was the same as previous studies (e.g., Gao et al., 2022;Song et al., 2022;Zheng et al., 2018).Water Resources Research 10.1029/2024WR037218

Validation of the Frozen Soil Module
Data sets related to frozen soil were collected to validate the model.The observed freezing and thawing front depth data at three China Meteorological Administration stations were collected (Zhiduo, Qumalai, and Yushu).The observed freezing/thawing front depths were compared with the simulated CF FS of the REW in which the station is located (Figure 12).The temporal variation of CF FS was consistent with the observed seasonality of frozen soil at all three stations.Specifically, the CF FS decreased during the freezing process (the frost depth increased) and increased during the melting process (the frost depth decreased and the thaw depth increased).The permafrost probability in each REW was extracted from the NIEER data set (Ran et al., 2022) to evaluate the spatial pattern of frozen soil simulation.As shown in Figure 13, the REW scale average CF FS had a significant negative correlation with the permafrost probability, indicating that the CF FS was lower in regions where more frozen soils existed, which was reasonable according to the definition of CF FS .However, it should be noted that the CF FS was just a coefficient introduced by the THREW-T model to correct the soil-related parameters based on the frozen soil conditions, and no equations related to the frozen soil characteristics were established.Consequently, these analyses could only provide a qualitative validation of the frozen soil module.

Discussions
The main limitation of this study was the simplification of the frozen soil module.First, the developed module just corrected two soil-related parameters according to the frozen soil conditions, but did not simulate the characteristics of frozen soil itself.The semi-distributed THREW-T model conceptualized the soil water and groundwater storage by two reservoirs (Tian et al., 2006), so it cannot simulate the detailed conditions of soil water and temperature at a soil profile like some other fully distributed models (e.g., Qi, Wang, et al., 2019;Song et al., 2020).Simulating the freezing and thawing processes was rather complex based on the structure of the current THREW-T model.Consequently, we developed a very simplified module, which only characterized the variation of parameters influenced by frozen soil.Second, when characterizing the hydrological impact of frozen soil, the simplified module only considered its influence on the soil properties (i.e., the soil-related parameters), but did not consider the influence on water amount (i.e., the state variables).In specific, the state variables related to soil water content (e.g., y s and y u in Equations 1 and 2) should be adjusted with the freezing and melting processes, and the melting water from frost could generate runoff directly (Gao et al., 2018;Song et al., 2022).However, the current module did not take into account the influence of water phase change when simulating the water balance.Finally, the CF FS was assumed to be corelated to the square root of accumulated ground temperature.Although the Equation developed in this study is similar to the form of Stefan equation, which is commonly adopted to calculate frost depth (Shur et al., 2005;Votyakova, 1976), there are many additional factors influencing the soil freezing and thawing processes, such as soil thermal conductivity, soil water content (Hayashi et al., 2007) and snowpack (Pomeroy et al., 2007;Qi, Wang, et al., 2019).The frozen soil module developed in this study adopted the spatially and temporally constant parameters to calculate CF FS , making the model unable to reflect these factors.Besides, due to the data scarcity on the TP, only data from two isotope sampling locations and three frozen soil measurement stations were available for model calibration and validation.The model representation over the whole region could be further improved by establishing more measurement sites.Although the proposed module was rather simplified, some indirect validation was conducted toward observation data sets, ensuring the reasonability of the frozen soil simulation.Nonetheless, further works are needed to incorporate modules with more physics into the model to improve its performance in cold regions, especially for the simulation of cryospheric processes.
Given the above limitations, the main aim of this study was to diagnose the hydrological model with tracer-aided method and improve the model performance in a cold basin, rather than developing a new frozen soil model.Diagnosing model structures and identifying the hydrological impact of some specific processes is an important way to improve the model performances (He et al., 2015;Tian et al., 2012), which is always difficult solely based on discharge data.For instance, although it is widely recognized that frozen soil has important influence on hydrological processes, many studies adopted models without frozen soil module in frozen soil influenced regions (e.g., Ahmed et al., 2022;Han et al., 2019;Mahmood et al., 2020).Even so, good discharge simulation could usually be produced by adjusting the parameters,  indicating that calibration solely toward discharge cannot fully capture the hydrological impact of some processes.However, our results demonstrated that by calibrating the model toward discharge and isotope data simultaneously, we were able to effectively identify the influence of frozen soil, ensuring that the model produced good results for right reasons.Thus, this study highlights the value of tracer-aided hydrological modeling method on diagnosing model structure by identifying the influence of specific processes such as frozen soil.

Conclusions
This study adopted a tracer-aided hydrological model THREW-T in a typical cold basin with widespread frozen soil on the TP.The physical representation of model structure was analyzed based on the shortage of simulation.The influence of frozen soil was identified and a module representing the hydrological impact of frozen soil was incorporated into the model.Our main findings are as follows: 1.The THREW-T model cannot perform well on the simulations of baseflow and stream water isotope simultaneously.The low discharge during dry seasons and the high contribution of subsurface runoff cannot be simultaneously characterized by the model.The imbalanced performance on two objectives could be attributed to the freeze-thaw cycle of frozen soil, resulting in seasonal variation of soil-related parameters, which was not considered in the THREW-T model.This highlights the value of tracer-aided modeling method on diagnosing model structure by identifying the hydrological impact of specific processes.2. A catchment-scale simplified frozen soil module was developed to correct the soil-related parameters (KKD and WM) according to the accumulated ground temperature.Incorporating the developed module into THREW-T model significantly improved the simulations on baseflow, isotope, and the balance of two objectives.The model incorporated with frozen soil module produced low baseflow and high contribution of subsurface runoff during wet seasons simultaneously.The introduced correction factor CF FS showed strong correlation with measurement front depth and permafrost possibility, validating the reasonability of frozen soil module.3. The frozen soil module revealed the hydrological effects of soil freeze-thaw cycle.The existence of frozen soil had little influence on the annual streamflow, but changed the runoff seasonality by reducing baseflow during dry seasons and increasing runoff during wet seasons.Meanwhile, the frozen soil also had little influence on the contribution of subsurface runoff in annual runoff, but increased this proportion during wet seasons.

Figure 1 .
Figure 1.Locations and topography the source regions of Yangtze River.

Figure 2 .
Figure 2. The spatiotemporal variation of precipitation δ 18 O in the source region of Yangtze River basin derived from the corrected isoGSM.(a) spatial distribution of average precipitation δ 18 O, (b) weighted average daily precipitation δ 18 O during the ablation period in 2016, (c) weighted average daily precipitation δ 18 O during the ablation period in 2017.

Figure 3 .
Figure 3. Schematic illustration of the THREW-T model.

Figure 4 .
Figure 4. Simulations of daily discharge, monthly snow cover area and daily stream water δ 18 O obtained by (a-c) NoFS-DS and (d-f) NoFS-DSI calibration variants.

Figure 6 .
Figure 6.The relationship between the representative elementary watershed scale mean annual ground temperature and air temperature.

Figure 7 .
Figure 7. Simulations of (a) daily discharge, (b) monthly snow cover area, and (c) daily stream water δ 18 O obtained by FS-DSI calibration variant.

Figure 8 .
Figure 8. Envelope line of the relation between the isotope and baseflow simulation obtained by NoFS-DSI and FS-DSI calibration variants.

Figure 9 .
Figure 9. Monthly discharge obtained by three calibration variants.(a) Monthly discharge of total runoff and subsurface runoff, (b) monthly contribution of subsurface runoff in total runoff (the solid and dashed lines represent total runoff and subsurface runoff, respectively).

Figure 10 .
Figure 10.Time series of basin scale average daily (a) ground temperature, (b) accumulated ground temperature, and (c) CF FS .

Figure 11 .
Figure 11.Comparison of simulated (a) daily discharge, (b) stream water isotope, (c) monthly discharge of total runoff and subsurface runoff, and (d) monthly contribution of subsurface runoff in NoFS and FS scenarios.

Figure 12 .
Figure 12.Time series of the daily simulated CF FS and observed freezing and thawing front depth at three stations.

Figure 13 .
Figure 13.The correlation between representative elementary watershed scale average CF FS and permafrost probability.

Table 1
Characteristics of Stream Water Samples and Isotope Composition in the Source Region of Yangtze River Basin

Table 2
Physical Descriptions of the Calibrated Parameters in the THREW-T Model NAN ET AL.

Table 3
Calibrated Parameter Values Obtained by Two Calibration Variants Figure 5.The daily isotope signature of multiple runoff components produced by the NoFS-DSI calibration variant.