Freezing Temperature Controls Winter Water Discharge for Cold Region Watershed

Pronounced climate warming over the arctic‐subarctic regions has lead to profound hydrological changes including intensified river flow, but how soil frost controls aquifer discharge remains poorly understood. This study quantifies the relationship between freezing temperature and baseflow in winter. Analyses show that the traditional reservoir models are unable to reproduce the observed baseflow variations. By incorporating a freezing temperature function in the reservoir models, the model performances are largely improved. It indicates the dominant role of freezing temperature in controlling the aquifer discharge through reducing the watershed conductivity and liquid (active) water content. The results for the Albany watershed in Canada show that the watershed lump conductivity decreases by half when air temperature accumulates to −172 °C·day from winter start and in extremely cold years, it could decrease by more than 85%. With this relationship, a climate warming of +1, +2, and +4 °C would suggest an increase of 7.7%, 16.7%, and 41.0% in conductivity or 6.8%, 14.7%, and 35.0% in winter discharge, respectively. The study provides an important link between climate warming and aquifer discharge in cold regions. The results could be particularly useful for developing process‐based models, estimating baseflow variations, and assessing climate change impact on cold region hydrology.


Introduction
Baseflow or recession flow is defined as the water depletion process of a watershed when there is no surface runoff or quick flow contributions to the river flows. Baseflow or recession analysis has been extensively used to understand many aspects of an aquifer, such as water storage, discharge mechanism, hydrogeological characteristics, and climate change impact, as reviewed in Hall (1968), Tallaksen (1995), and Smakhtin (2001). With the terrestrial total water storage (TWS) made available by the GRACE satellites, the baseflow analysis has been extended for estimating drainable water storage and quantifying storage-discharge relationships (Riegger & Tourian, 2014;Sproles et al., 2015;Tourian et al., 2018;Wang et al., 2017;Wang & Russell, 2016).
Recession models have been commonly derived from the basic differential equation governing discharge from an unconfined aquifer. Brutsaert and Nieber (1977) demonstrated that the general recession behavior can be studied by characterizing the baseflow with the power-law model: where Q(t) is the baseflow of an aquifer at time t and c and d are constants. The values of c and d can be obtained by plotting log(−dQ/dt) against log(Q), where the slope is d and the intercept is log(c). Typically, the time increment dt over which the recession slope dQ/dt is calculated is held constant. However, Rupp and Selker (2006) showed that the use of constant dt could limit accurate representation of the recession relationship to the portion of the hydrograph for which the chosen dt is inappropriate. It could also lead to artifact upper and lower envelopes in the graphs. Biswal and Marani (2010) and Kumar (2014a, 2014b) proposed an alternative geomorphological recession flow model to improve recession analysis for sloped watersheds. Their model considers the temporal evolution of active drainage network length for recession flows and has demonstrated significant improvement in estimating the recession curve properties.
The above power-law recession model is consistent with a class of reservoir models, which can be generalized as where S(t) is aquifer water storage at time t and k and b are parameters representing the watershed lump hydraulic conductivity and storage-discharge relationship, respectively. When b = 1, the storage-discharge has a linear relationship. Otherwise, the model becomes nonlinear (b < 1 less linear; b > 1 more linear).
In equation (2), k = [c(2 − d)] 1/(2 − d) , and b = 1/(2 − d) (Clark et al., 2009;Kirchner, 2009), such that the linear recession model d = 1 corresponds with the linear reservoir model b = 1. Available studies found that b is rarely less than 1 and values greater than 1 are typical for many watersheds (Stoelzle et al., 2013;Tallaksen, 1995). While the assumption that b = 1 is not applicable for all applications, Vogel and Kroll (1992) found this to be a valid assumption in selected northeastern United States watersheds. The study by van Dijk (2010) found in 183 Australian catchments that a linear reservoir model produced baseflow estimates as good as those obtained using a nonlinear reservoir and for practical purposes a linear reservoir model was therefore recommended.
Parameter k is commonly regarded as being determined by watershed physical and geological characteristics, including geomorphology (e.g., watershed size, slope, and drainage density) and soil and aquifer characteristics (e.g., hydraulic conductivity, porosity, or specific yield; Farvolden, 1963;Bloomfield et al., 2009;Stoelzle et al., 2014;Chiverton et al., 2015;Paznekas & Hayashi, 2016). Brandes et al. (2005) found that geomorphological indices such as drainage density, slope, and hydrologic soil class together explained about 70-80% of the variation in k for 24 watersheds in the Appalachians (USA). This makes it possible to derive average values of k for different physiographical regions or geological formations and apply them as indices of representativeness in hydrological regions. The hydroclimatic conditions and the human interference in the water cycle (e.g., groundwater withdrawal) during the recession process are also found to affect the estimation of k (Bogaart et al., 2016;Patnaik et al., 2018;Price, 2011). For example, van Dijk (2010) related k to climate factors including precipitation and evapotranspiration and found that 27% of the variance in k could be attributed to variability of humidity for arid watersheds. Thomas et al. (2013) described the impact of human interference of groundwater withdrawal on the estimation of k. The values of parameter b also vary with watershed characteristics and hydroclimatic conditions. Tallaksen (1995) suggests that a continually changing relationship between storage and baseflow can be accounted for by the changes of b values. These findings provided insights on the impact of water conditions on watershed discharge modeling.
Baseflow in cold regions is more complicated due to the seasonal cycles of soil freeze and thaw. Soil frost affects the baseflow through a number of processes such as limiting water infiltration through the vadose zone for aquifer recharge, reducing the liquid water storage capable for discharge, altering watershed hydraulic conductivity, and even changing the local water flow paths in a watershed (Ireson et al., 2013). Even in summer, relatively cold antecedent winter soils were found to have a propagating effect on the recession (Ploum et al., 2019), suggesting that frozen soil layers present in the vertical soil profile can affect recharge of deep groundwater. Baseflow analysis traditionally focuses on linking recession parameters k and b with aquifer physical or geological characteristics that are often thought to be relatively constant (Seibert & Meerveld, 2016;Troch et al., 2013). Since seasonal soil freezing is a temporarily dynamic phenomenon determined mainly by external factors such as air temperature, the impact of soil frost on discharge parameters k and b has rarely been studied and remains poorly understood.
Cold regions in high latitudes have experienced stronger warming than lower latitudes in the past decades (Bekryaev et al., 2010;Francis et al., 2017). The pronounced warming has been found to alter many aspects of the cold region hydrological process such as streamflow magnitude and seasonality, partitioning of water into surface and subsurface, dormant aquifer activation, and even disappearing of lakes due to permafrost thaw (Evans et al., 2018;Green et al., 2011;Kurylyk et al., 2014;Lyon & Destouni, 2010;Smith et al., 2005;Taylor et al., 2013;Walvoord & Kurylyk, 2016). Durocher et al. (2019) analyzed a total of 72 rivers that contribute flow to the Arctic Ocean and discovered a general increasing trend of flow in the past four decades. Understanding and quantifying the relationship between temperature and the discharge rate in these regions are important for assessing and projecting climate change impact on baseflow. Since baseflow water temperatures from underlying aquifers may be different from the surface water such as snowmelt in spring (Hayashi & Rosenberry, 2001), the change in baseflow may affect river ice dynamics such as the ice formation and breakup. Spring flood from snowmelt is a common hazard for many rivers in cold regions, and it is mainly a result of ice jam (Corston, 2016;Rokaya et al., 2018). Traditionally, river ice dynamics has been perceived as being controlled by surface water. There is however growing recognition that winter baseflow can have significant implications for the river ice dynamics and spring flood.
The main objective of this study is to quantify the impact of freezing temperature on baseflow and to examine the storage-discharge relationship so that baseflow estimates can be improved for cold region watersheds. Although Moore et al. (2002) and Thomas et al. (2013) have noticed the recession rate variations with temperature, this study is the first one to quantify the dynamic changes of the baseflow parameters (k and b) with the freezing temperature in winter. New storage-discharge models are proposed to incorporate soil frost impact on aquifer discharge based on accumulated temperature in winter. This study will also examine the linearity (or nonlinearity) of the storage-discharge relationship and how this relationship will change with the freezing temperature. The rest of this paper is organized as follows: section 2 describes the hypotheses, model experiment design, and test methods; section 3 describes the study region and data used in this study, including the uncertainties analyses; the results are provided in section 4 and discussed in section 5, followed by conclusion remarks in section 6.

Hypotheses and Methods
Based on the above discussions on soil frost impact on the cold region hydrological processes, this study hypothesizes that freezing temperature can change aquifer discharge through changing the parameter values of k or b in the baseflow model (equation (2)). This study adopts the effective water storage concept as proposed in Wang et al. (2017) in the baseflow modeling. With these hypotheses and modifications, equation (2) can be written as where Q(t) in this study is expressed in depth of water per unit time (length time −1 ), S(t) is in depth of water (length), k 0 (length (1 − b) time −1 ) and b 0 (dimensionless) are the respective base values of parameters k and b when soil frost is not present, f(T) (dimensionless) is a function of air temperature T (°C), and a (length, water depth) is a parameter representing the threshold value of water storage below which the aquifer discharge would be 0. Equation (3a) sets k dynamic with f(T), hypothesizing that soil frost affects discharge through reducing the hydraulic conductivity of an aquifer. Equation (3b) sets b dynamic with f(T), hypothesizing that soil frost affects discharge through reducing the liquid (active) water content of an aquifer. Since soil frost is a result of both the current and antecedent air temperatures, accumulated temperature is used to formulate f(T): where T day is daily air temperature (°C), T acc (t) (°C·day) is the accumulated T day from the starting date of winter (t 0 ) to day t, t end is the ending date of winter, and T c (°C·day) is a parameter. The determination for t 0 and t end is given below.
The T acc represents the accumulated coldness at a specific date since the winter start. It is used as a proxy for soil freezing condition and has the advantage of easy to obtain. To my knowledge, this is the first time to use T acc for characterizing soil frost conditions. Equation (4a) is similar to the Michaelis-Menten function (Michaelis & Menten, 1913), which reduces k 0 by half when T acc reaches the value of T c . The function is used to describe the relationship between soil freezing conditions and the accumulated coldness at a specific date. An example for this relationship simulated by the EALCO model (Wang et al., 2013(Wang et al., , 2015Zhang et al., 2008) is shown in Figure S1 in the supporting information.
The above models contain up to four unknown parameters (k, T c , b, and a), depending on different hypotheses. The values for these parameters were found by solving a nested numerical iteration scheme with the criteria of maximum Nash-Sutcliffe modeling efficiency (NSE) for baseflow. Model experiments with five different hypotheses were conducted. Table 1 summarizes the model hypotheses and parameters to be solved.
This study is only concerned with the winter season, during which the river flow is mainly contributed by aquifer discharge or baseflow. The start of the winter season (t 0 ) is determined by the criteria of (1) the daily average temperature drops below 0°C in the late fall and (2) the accumulated temperature after this date remains negative. The second criterion is to exclude temporary cold spans that would not lead to the soil in a continuous freezing conditions in the early winter season. The end of the winter season (t end ) is determined by the criteria of (1) the daily average temperature rises above 0°C in the spring and (2) the accumulated temperature after this date remains above 0. The second criterion is applied to exclude the possible cases that the temperature momentarily rises above 0°C but would not result in prolonged land surface thaw. The determination of the "winter season" is quite conservative in order to avoid time periods with possible contributions from surface runoff (e.g., rain or snowmelt) to the observed flows.
The models were tested using monthly baseflow values due to the facts that (1) the GRACE temporal resolution is monthly and (2) since the watershed is large and the major river is about 1,000 km long, it takes time for the discharged water to flow to the gauge station at the river mouth. So there would be a time lag between discharge and gauge measurement. Also, there would be a smoothing process when the water flows down to the river mouth. As such, it would be problematic to compare modeled discharge and gauge station measurement at a daily time step. Monthly time step reduces the impact of the above-mentioned processes on baseflow modeling, although it does not completely remove the problems. The water storage at month m (S m , mm) was obtained from its value in the previous month S m − 1 (mm) and the amount of water loss from observed baseflow in that month. The loss of water from soil and aquifers due to evapotranspiration is ignored due to the snow-covered ground. The water storage at the beginning of a winter season S 0 was obtained from the GRACE TWS data set as detailed in the next section.
The model performance was evaluated using the mean absolute error (MAE), the Pearson correlation coefficient (R) and its significance test (p), and the Nash-Sutcliffe model efficiency coefficient (NSE, equation (5)).
The NSE is a normalized statistic that determines the relative magnitude of the residue variance compared to the measured data variance. It is commonly used to assess the predictive power of hydrological models. Model Note. Five models (Exp-1 to Exp-5) with different hypothesis are tested (Q: baseflow; T: air temperature; S: water storage; definitions and dimensions for the parameters in the table can be found in section 2).

Study Region and Data Sets
The above hypotheses were tested for the Albany watershed located in the cold region of far north of Ontario, Canada ( Figure 1). It was selected based on the facts that it is the largest watershed in Ontario (with a drainage area of 137,230 km 2 ) to fit with the large footprint of GRACE satellites, has gauge measurement available for the study period, and is highly sensitive to climate change (McLaughlin & Webster, 2014). The watershed is highly vulnerable to flooding especially during snowmelt season in spring due to ice jam. The Albany River has a length of 982 km, flowing northeast from Lake St. Joseph at an elevation of 371 m in Northwestern Ontario into James Bay. The headwater of the river basin is situated in the Canadian Shield physiographic region, which is characterized by a thin soil layer over Precambrian bedrock and moderate topographic relief. The middle and lower portions of the drainage basin are within the Hudson Bay Lowlands physiographic region, which is characterized by poorly drained organic deposits and low topographic relief. Within the study area, the Canadian Shield landscape is dominated by Boreal forest, which transitions into Barren Boreal and Taiga vegetation zones within the Hudson Bay Lowlands.
Three data sets were used, which include air temperature, river flow, and GRACE TWS. Limited mainly by the record length of GRACE satellites data, the hypotheses were tested over the time period of 2002-2015, which includes a total of 13 winter seasons. These data sets and the regional hydroclimatic conditions are described below.

Air Temperature (T)
The daily T for the watershed was calculated from the Global Land Data Assimilation System (GLDAS) meteorological forcing. The GLDAS data were acquired as part of the mission of NASA's Earth Science Division and archived and distributed by the Goddard Earth Sciences Data and Information Services Center (http://mirador.gsfc.nasa.gov/). The data are provided at a 3-hr time step and a spatial resolution of 0.25°by 0.25°latitude/longitude. The daily T of a grid is taken as the average of the eight readings of the day, and the watershed T is the average of the T values of all the GLDAS grids within the watershed. There are four weather stations located in the Albany watershed with air temperature measurement available during the study period. A comparison of the station measurements and the GLDAS temperatures in the corresponding grids shows that they are highly consistent, with average biases varying from −0.85 to 0.01°C and linear correlation coefficients varying from 0.96 to~0.99 ( Figure S2).
Study area temperature has variations from cold winters to mild summers. In the study period, the temperature on average dropped below 0°C in early November and rose above 0°C in middle April of next year. Both of the transitional time were of large inter-annual variations for more than a month. Air temperature mostly reached the lowest of −25 to −30°C in late January. The accumulated temperature T acc in the winter season varied between −1,250 and −2,080°C·day ( Figure 2). The extremely low temperature and low solar radiation in the winter season result in the deep frozen condition for the watershed in winter.

River Flow Measurement (Q)
The Albany River has several hydrometric stations managed by Water Survey of Canada. This study used data from the most downstream gauge station in order to have the maximum drainage area possible (Gauge Station #04HA001, Albany River near Hat Island, latitude 51.3306, longitude −83.8333). The observed daily Q at the hydrometric station was downloaded from the Water Survey of Canada website (http://www.ec.gc.ca/rhc-wsc/). The original Q is in m 3 s −1 , and it was converted to water depth (mm) using the drainage area. The quality of the Q data was cross checked by inter-comparisons with nearby gauge measurements over the region. Specifically, there are two major watersheds adjacent to the Albany watershed: the Attawapiskat watershed in the northwest (Station #04LG004) and the Moose watershed in the southeast (Station #04FC001; Figure 1). Comparisons showed that the Q for Albany was reasonably correlated with those for Attawapiskat and Moose, with a linear correlation coefficient of 0.84 and 0.78 ( Figure S3), respectively. No suspicious outliners were found in the data sets. The Q data were also compared with two gauge measurements within the Albany watershed, which represent the two major contributory sub-watersheds, namely, the Kenogami River (Station #04JG001) and the Albany River above Nottik Island (Station #04GD001). The linear correlation coefficients with these two sub-watersheds were 0.91 and 0.83 ( Figure S3), respectively. The above watersheds are located in a large flat region with similar synoptic and land surface conditions. The correlations among these gauge measurements are expected, and they demonstrate the consistency and reliability in the Q data quality.
The hydrology of the Albany River is characterized by an annual mean discharge of 1,420 m 3 /s (380 mm year −1 ) with peak flows as high as 8,000 m 3 /s (5.86 mm day −1 ) induced by snowmelt in April to late May.
Watershed discharge decreases continually in the winter season, and it presents a typical recession process since there is little contribution from surface runoff to the river flow (Figure 3). The inter-annual differences in the baseflow values are large in early winter, varied from above 4,600 m 3 /s (3.37 mm day −1 ) to just 350 m 3 /s (0.26 mm day −1 ) during the study period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). The minimum baseflow values in late winter . It is worth mentioning that the annual precipitation amount for the Albany basin is around 700 mm and the annual evapotranspiration is less than half of it (Wang et al., 2013), so the basin has a large water surplus for aquifer recharge to sustain the winter baseflow.
Plot of log(−dQ/dt) versus log(Q) for the Albany watershed in the winter months of 2002-2015 was given in Figure 4. First, a large intra-season change in the slope (d) can be found in each year. In early winter, the slope remained relatively constant. Along with the progress of the winter season, the slope increased significantly, resulting in the large decrease in the value of intercept (log(c)). This large intra-season change is consistent for each year, suggesting that the baseflow model parameters are dynamic and affected by factors

Watershed Water Storage (S)
The watershed water storage at the beginning of winter S 0 was estimated using the TWS product from GRACE Release-05 spherical harmonic solutions. The data were downloaded from the GRACE Tellus website (ftp://podaac-ftp.jpl.nasa.gov/allData/tellus/L3/land_mass/RL05/), which is supported by the NASA MEaSUREs Program (Swenson, 2012). The GRACE models contain recoverable gravity change signals up to a maximum spherical harmonic degree 60, which approximately represents a spatial resolution of~330 km (108,900 km 2 ). The land grid scale factor, a set of scaling coefficients intended to restore much of the energy removed by the destriping, Gaussian, and degree 60 filters to the land grids, was applied to the data sets  to minimize the difference between the model's smoothed and unfiltered monthly water storage variations at any geographic location (Landerer & Swenson, 2012). The data sets include monthly TWS time series from three data processing centers: namely, University of Texas/Center for Space Research (CSR), GeoForschungsZentrum Potsdam (GFZ), and Jet Propulsion Laboratory (JPL). The differences of the three TWS data sets over Canada's landmass were investigated in Wang and Li (2016). The differences were found minor for the Albany watershed (average absolute difference was 9 mm). In this study, the average TWS of the three data sets was used. The baseline of the TWS, which was based on the average from January 2004 through December 2009 in the original data, was readjusted to the minimum value found over the study period.
The TWS for the Albany watershed had an average seasonal variation range of 100 mm during the 13-year study period ( Figure 5). The lowest values appeared in late summer around September and October, and the highest values occurred in the late spring around April before snowmelt began. Obviously, the increase in TWS in the fall-winter season is mainly due to the snow accumulation and low evapotranspiration, and the decrease in the spring-summer season is mainly due to the large amount of discharge of snowmelt water and the high evapotranspiration in summer. The inter-annual variations in TWS for a specific month were large, varied from 70 mm appeared in June to 127 mm appeared in November. Monthly measurement error, leakage error, and combined total error in the GRACE TWS were calculated following Wahr et al. (2006), and they were 13.2 mm, 15.8 m, and 20.6 mm, respectively, for the watershed. The impact on the TWS error estimate due to the uncertainty in the CLM4 land surface model used for the scale factor calculation was evaluated by comparing with a different land surface model of NOAH in an independent study (Wang et al., 2014). The magnitudes of errors from the two land surface models were found to be similar.

Results
The traditional linear reservoir model without having the temperature function (Exp-1) was unable to reproduce the observed recession process. Specifically, the model explained only 12.2% of the variations in the observed monthly baseflow values in the 13 winter seasons ( Figure 6). The model performance test showed a low NSE of 0.122 (Table 2).
The traditional nonlinear reservoir model without incorporating the temperature function (Exp-2) achieved the best model performance with the exponent factor b = 1.3 (Table 2), suggesting a nonlinear storagedischarge relationship. However, the nonlinear model only slightly improved the linear model performance in Exp-1. Specifically, the nonlinear model explained 13.1% of the variations in the observed baseflow values ( Figure 6). It is only slightly higher than that obtained in the linear model (12.2%). The limited improvement is also reflected by the small increases in the NSE ( Table 2). The poor model performance in both Exp-1 and  Table 1.
Exp-2 suggests that water storage is not the dominant factor controlling the discharge process for the watershed in winter.
With the modified linear reservoir model (Exp-3), which incorporated the freezing temperature function in the hydraulic conductivity, the model performance was largely improved as can be seen from Figure 6. The model explained 81.6% of the variations in the observed baseflow values in the 13 winter seasons. The R between the modeled and observed monthly baseflow reached above 0.9 (p < 0.001). The mean absolute errors (MAE) were decreased by more than a half compared with that in Exp-1 and Exp-2. The improvement can also be seen from the high NSE value of 0.809. This large improvement in model performance clearly demonstrates the dominant role of freezing temperature in the discharge process. The parameter value for T c in Exp-3 indicates that the k was decreased by 50% when the accumulated temperature reached to −171°C·day from the winter start.
By incorporating the temperature function in the hydraulic conductivity in the nonlinear reservoir model (Exp-4), the results were further improved over that of Exp-3, but the improvement is small (Table 2). This is consistent with the comparison between Exp-1 and Exp-2, which suggests that employing nonlinear storage-discharge relationship, either with static k (Exp-1 and Exp-2) or dynamic k (Exp-3 and Exp-4), has limited significance in improving the baseflow modeling in winter for the Albany watershed. The above results are based on static b. It is interesting to note that the parameter b has a value of 1.01 in Exp-4, which is not significantly different from the b in Exp-3 (=1.0) when assuming an uncertainty level of 10% in flow measurements (Wang et al., 2014). As such, Exp-4 virtually converged with Exp3, and it suggests a linear storage-discharge relationship when k is set dynamic. The parameter T c in Exp-4 has a value of −172°C·day, which is consistent with that obtained in model Exp-3. For model simplicity, the result from Exp-3 is used in the illustrations hereafter.
The dynamic k values obtained in model Exp-3 and inversely calculated using the observed baseflow for all the months in the study period have a Pearson correlation coefficient close to 0.9 (p < 0.001; Figure 7). The high correlation demonstrates the dominant role of freezing temperature in controlling the variations of hydraulic conductivity for the watershed in winter. The T-k relationship determined by model Exp-3 is compared well with that inversely calculated using the observed baseflow ( Figure 8). It shows that the k was decreased from its pre-winter value of 0.15 month −1 to as low as 0.022 month −1 , or by about 85%, in extremely cold winters when the accumulated temperature reached to about −2,000°C·day late in the season. With this relationship, a climate warming of +1, +2, and +4°C would suggest an increase of 7.7%, 16.7%, and Note. Exp-1: traditional linear reservoir; Exp-2: traditional nonlinear reservoir; Exp-3: linear reservoir with dynamic k; Exp-4: nonlinear reservoir with dynamic k; and Exp-5: nonlinear reservoir with dynamic b. The detailed modeling hypotheses for Exp-1 through Exp-5 can be found in Table 1 (Mean = mean baseflow, MAE = mean absolute error, R = Pearson correlation coefficient, NSE = Nash-Sutcliffe model efficiency, p = significance level).

Figure 7.
Comparisons of the baseflow parameter k estimated using observations (kobserve) and the model Exp-3 (kmodel). Exp-3 is a linear reservoir model with dynamic k (see details in Table 1). 41.0% in the k values, respectively. This would result in an increase of 6.8%, 14.7%, and 35.0% in winter water discharge or 3.53 × 10 8 , 7.60 × 10 8 , and 1.81 × 10 9 m 3 by water volume, respectively, for the three warming scenarios for the Albany watershed. It is worth noting that these estimates are based on the impact of temperature on discharge alone, which assume that other climate conditions of the watershed remain the same.
When the exponent factor b was set to dynamic with the freezing temperature function (Exp-5), the model performance also showed large improvement over the traditional linear and nonlinear reservoir models of Exp-1 and Exp-2 ( Figure 6). Specifically, the results from Exp-5 explained 79.6% of the variations in the observed baseflow values in the 13 winter seasons. The NSE reached a value of 0.782. The results demonstrate that the working hypotheses of soil frost affecting aquifer discharge through reducing the liquid (active) water content are valid. It is worth noting that that the b 0 value obtained in Exp-5 is 1.0, suggesting a linear storage-discharge relationship if soil frost is absent (i.e., f(T) = 1.0) in the watershed. This is consistent with the findings from Exp-3 and Exp-4. The Exp-5 slightly underperformed Exp-3 and Exp-4. However, the largest difference in the model performance occurred between the models with static k and b (i.e., Exp-1 and Exp-2) and the models with dynamic k and b (i.e., Exp-3, Exp-4, and Exp-5). This clearly demonstrates that the parameters k and b in the baseflow model (equation (2)) are actually dynamic in nature and they are strongly affected by the soil freezing conditions in winter.
The dynamic b values obtained in model Exp-5 and inversely calculated using the observed baseflow for all the months in the study period have a Pearson correlation coefficient of 0.876 (p < 0.001; Figure 9). The T-b relationship obtained by model Exp-5 and inversely calculated using the observed baseflow is compared in Figure 10. The range of b values during the study period varied from its pre-winter value of 1.0 when there is no soil frost to a lowest value of 0.53 in the late winter of 2013 and 2014 when the T acc reached to the lowest value of −2,000°C·day.

Discussion
Baseflow analysis has been facing numerous challenges. First, it is important for a model to include the dominant factors controlling the discharge process. In this study, the traditional reservoir models (Exp-1 and Exp-2) were found unable to reproduce the baseflow observed in winter in the Albany watershed, suggesting that the aquifer discharge in cold regions was mainly controlled by physical factors apart from water storage. By using a freezing temperature function to modify the reservoir models (equation (2)), large model improvement was achieved. It demonstrates that freezing temperature is a dominant factor in controlling aquifer discharge in winter. The controlling process can be modeled through setting either k or b dynamic with the temperature function or through reducing either the hydraulic conductivity or the liquid (active) water content with freezing temperatures. Another important finding is that by incorporating the freezing temperature function in the reservoir models, the watershed storage-discharge relationships in Exp-3, Exp-4, and Exp-5 were converged to b 0 = 1.0, suggesting a linear relationship when soil frost is absent. The linear storage-discharge relationship obtained for the Albany watershed is not surprising. The watershed is located in a flat region with a large water budget surplus for aquifer recharge before the winter season. The winter baseflow is relatively simple and mainly  Table 1).

Figure 9.
Comparisons of the dynamic values of the exponential factor (b) in the baseflow model (equation (2)). bobserve: estimated using observations, bmodel: estimated from the model Exp-5. Exp-5 is a nonlinear reservoir model with dynamic b (see details in Table 1).
contributed by surficial aquifers including wetlands. The hydrogeological settings of the watershed are largely different with the watersheds that have deep slops and complicated aquifer systems with large variations in water conditions, which often involve dynamic channel networks (e.g., time-varying geometry of saturated channeled sites) and result in nonlinear storage-baseflow relationships (Biswal & Kumar, 2014b;Biswal & Marani, 2010;Tallaksen, 1995).
Second, baseflow analysis is often challenged by the selection of characteristic recession segments. The determination of the recession segments imposes large variabilities in the later modeling of baseflow. This is because the selected recession segments may represent different water conditions prior to the recession, different stages in the temporal evolution of its active drainage network, and different climate conditions during the recession process. For example, numerous studies have reported that the estimated recession parameters depend on both the start and length of the recession segment and the seasonal change in evapotranspiration (e.g., Biswal & Kumar, 2014a;Federer, 1973;Ficklin et al., 2016). During warm and active plant growing season particularly for regions with shallow groundwater tables where a drying of the upper soil layer can be succeeded by capillary transport from the groundwater, evapotranspiration could result in a reduction in baseflow and overestimate in c in equation (1) or k in equation (2) (Federer, 1973;Tallaksen, 1995). In many cases, the selected recession segments involve water contributions from surface runoff, so a baseflow-surface runoff separation from the total flow measurement is required. This often brings another source of large uncertainties in the baseflow data. This study overcame the above-mentioned challenges by selecting winter season for the recession analyses. First, surface runoff is minimal in winter due to the lack of liquid precipitation and snow cover, so the winter flow was essentially contributed by aquifer discharge. It avoids the baseflow-surface runoff separation process and therefore the large uncertainties associated with it. Second, the recession segments selection using air temperature is rather objective and less dependent on users than the selection process from hydrography records. Third, there is minimal impact by evapotranspiration on the recession process as in the winter season the ground surface is covered by snow and groundwater loss through evapotranspiration is minimal. The aquifer discharge process in winter is also largely simplified by the reduction of exchange between surface water and groundwater due to the frozen soil. Moreover, the long winter season for cold region watersheds provides the advantage of maximal data availability for model test. This study uses GRACE TWS for initial water storage estimation. Determining large-scale aquifer water storage is extremely difficult by traditional methods. This study demonstrates an innovative application of GRACE satellite observations in aquifer discharge and recession studies.
It is well known that different baseflow models often result in considerable variations in the estimates for model parameters or even different storage-discharge relationships for the same watershed (e.g., linear vs. nonlinear; Patnaik et al., 2018;Stoelzle et al., 2013;Tallaksen, 1995). This can be clearly seen by comparing the results from the five model experiments. For example, the k values varied from a low of 0.013 month −1 in Exp-2 to a high of 0.276 month −1 in Exp-3 (Table 2). However, as mentioned in section 1, the k is traditionally regarded as a physical parameter that represents the watershed, soil, and aquifer characteristics, which are often thought to be relatively constant. The large variations of k among the different models clearly indicate that it is actually dynamic in nature. For the cold region watershed in this study, the dynamics of k is found to be strongly related to the soil freezing conditions in winter, which decreased from its pre-winter value when soil frost was absent by about 85% in late winter when the soil was in deep frozen condition. Without including the dynamic variations of freezing temperature in the reservoir model (Exp-1 and Exp-2), the static k values represented a lump average of the dynamic k values during the entire winter season. As such, they were much smaller than the k 0 values obtained in Exp-3, Exp-4, and Exp-5, which represented the pre-winter condition when there was no soil frost. This further illustrates the importance for a model to include the dominant factors controlling the discharge process.  (2)) with accumulated temperature in winter (dots: estimated using observed baseflow, line: estimated using model Exp-5). Exp-5 is a nonlinear reservoir model with dynamic b (see details in Table 1).

Water Resources Research
WANG 10,490 Baseflow models can be tested using either a single recession segment or a group of recession segments from different time and conditions. Biswal & Marani, 2010, Biswal & Kumar, 2014b demonstrated that in the case when dQ/dt is dominated by the change of the length of active drainage network, model analysis based on a single recession segment can avoid a severe underestimation of the exponent factor from binning all recession segments together. Also, model test based on single recession segment is often easy to achieve relatively good results. Indeed, when the traditional reservoir models were applied separately to each of the winter seasons in this study, their performances were mostly better than those of Exp-1 and Exp-2. However, large variations in the values for the baseflow model parameters among the different years were noticed, which is not surprising due to the variations in climate conditions for different years. On the contrary, model test based collectively on all the recession segments imposes great challenge due to the impact of various factors on the recession process as discussed above. This challenge is also an opportunity to help discover important physical factors or mechanisms controlling the variations in the observed discharge values. In this study, the model tests using collectively all the data from the 13 winter seasons helped find the problems with the traditional reservoir models and reveal the dynamic impact of freezing temperature on the discharge process. The temperature function used in the baseflow model effectively consolidated the large variations in the model parameter values among different years. It also resulted in the convergence to a linear reservoir model for the Albany watershed when soil frost is absent. The findings demonstrate the importance of modeling approaches in explaining the variabilities of the observations from varying conditions and time. The findings can also help make more meaningful comparisons across different studies for different regions. The results are quite encouraging given the performance in the model test and the length of study period. As found in Perzyna (1993), the estimates of the recession parameters can be deemed as reliable when the study period reaches about 10 years.

Conclusions
This study discovered and quantified the impact of freezing temperature on winter discharge for a cold region watershed in Canada. It provides an important link between climate warming and baseflow change in cold regions. The results show that the traditional linear and nonlinear reservoir models were unable to reproduce the observed baseflow variations in winter, and they suggested that water storage is not the dominant factor controlling the discharge process. By incorporating a freezing temperature function in the reservoir models, the model performances were largely improved. This large improvement demonstrates the dominant role of freezing temperature in controlling the aquifer water discharge. This controlling process can be modeled through reducing the hydraulic conductivity of the aquifer or, specifically, through reducing the values of slope parameter k in the reservoir models. The results indicate that the k was reduced by half when the temperature accumulated to −172°C·day after winter start for the Albany watershed. In some extremely cold winters during the study period, the k was decreased by 85% when the accumulated temperature reached to −2,000°C·day in the late season. The controlling process of freezing temperature in aquifer discharge can be also modeled through reducing the liquid (active) water content or, specifically, through reducing the values of the exponent parameter b in the reservoir models. With this hypothesis, the model experiment shows that the b values varied in a range from its pre-winter value of 1.0 to as low as 0.53 in late winter during the study period. Another important finding is that by setting the parameters (k or b) in the reservoir models dynamic with the freezing temperature function, the watershed storage-discharge relationships were converged to a pre-winter base value of b 0 = 1.0, suggesting a linear relationship when soil frost is absent. With the temperature relationship found in this study, a climate warming of +1, +2, and +4°C would suggest an increase of 7.7%, 16.7%, and 41.0% in hydraulic conductivity or an increase of 6.8%, 14.7%, and 35.0% in winter water discharge, respectively. This is equivalent to a volume of 3.53 × 10 8 , 7.60 × 10 8 , and 1.81 × 10 9 m 3 more freshwater input into the James Bay from the Albany River, respectively. This study used the Albany watershed for demonstration based on considerations of data suitability and availability, sensitivity to climate change, and importance in social economics, but the hypothesis and approach are not watershed specific. In applications, the approach could be constrained by GRACE data resolutions when a watershed size is small. Another limitation is that the GRACE data are available only from 2002, which makes the study period relatively short. With the continuation of the GRACE Follow-On mission and the advancement in data processing, the above limitations are expected to be improved. Nevertheless, this study provides insight into the freezing temperature impact on cold region aquifer