Implementation of salt-induced freezing point depression function into 1 CoupModel _ v 5 for improvement of modelling seasonally frozen soils 2

13 Soil freezing/thawing is important for soil hydrology and water management in cold regions. Salt in 14 agricultural field impacts soil freezing/thawing characteristics and therefore soil hydrologic process. In 15 this context, we conducted field experiments on soil water, heat and salt dynamics in two seasonally 16 frozen agricultural regions of northern China to understand influences of salt on cold regions hydrology. 17 We developed CoupModel by implementing impacts of salt on freezing point depression. We employed a 18 Monte-Carlo sampling method to calibrate the new model with field observations. The new model 19 improved soil temperature mean error (ME) by 16% to 77% when new freezing point equations were 20 implemented into CoupModel. Nevertheless, we found that parameters related to energy balance and soil 21 freezing characteristics in the new model were sensitive to soil heat and water transport at both sites. 22 * Correspondence author. Email: jingwei.wu@whu.edu.cn Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess-2018-466 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 8 October 2018 c © Author(s) 2018. CC BY 4.0 License.


Introduction
Soil freezing and thawing processes have long been recognized for its importance in not only engineering applications (e.g., construction of roads and pipelines) (Jones, 1981;Hansson et al., 2004;Wettlaufer and Worster, 2006), but also environmental issues (e.g., soil erosion, flooding, and pollutants migration) (Andersland et al., 1996;Seyfried and Murdock, 1997;Baker and Spaans, 1997 ;McCauley et al., 2002).Knowledge on soil freezing and thawing could uncover mechanisms on water and salt distribution in soil (Baker and Osterkamp, 1989), on frost heaving (Wettlaufer and Worster, 2006), on waste disposal technology (McCauley et al., 2002), as well as on climate change and water management in cold regions (Lopez et al., 2007).
Laboratory and field experiments on hydrological characteristics of freezing/thawing soils have been conducted to understand soil hydrology in cold regions.Most of the experiments focused on soil freezing characteristics under various climate and soil conditions (Williams, 1964;Black and Tice, 1989;Spaans and Baker, 1996;Azmatch et al., 2012), regional water and energy balance in winter (Fuchs et al., 1978;Baker and Spaans, 1997;Hayashi et al., 2004;Watanabe et al., 2013;Zhou et al., 2014).There were very few studies on salt transport in frozen soils, except for frost heaving.Cary et al. (1979) found salt can decrease frost heaving and increase infiltration in frozen soils based on observations.Konrad and McCammon (1990) found the expulsion of salt from ice is dependent on freezing rate of soil.Hydrological effects of salt in cold regions have not been deeply explored.Wang et al. (2016) compared water and salt fluxes in two agricultural fields same as in this study, and detected different flow characteristics of salt during soil freezing and thawing seasons.They demonstrated that salt expulsion and dispersion are not negligible in frozen soils.Wu et al. (2016a) found that evaporation during winter was controlled by soil salt and groundwater in field frost tube experiments, in Inner Mongolia, China.They also demonstrated that water, heat and salt transport in frozen soils were coupled, and due to spatial heterogeneity of soil properties and technical difficulties in soil freezing/thawing experiments, measurements contained large uncertainties.
Numerical models on soil freezing/thawing have been put forward by many.Jansson and Karlberg (2004) developed a coupled process-based model-CoupModel, to simulate water, heat as well as salt transport in frozen soil.CoupModel is a process-based model with detailed descriptions on coupled water and heat transport in frozen and unfrozen soils.It has shown to be one of the most robust models among other models taking soil freezing/thawing into account (e.g., SWAP, DRAINMOD, SWAT, HBV, VIC, and ATS etc).This model was developed and applied to forests (Gustafsson et al., 2004;Wu et al., 2013), agricultural field (Wu et al., 2011), permafrost (Zhang et al., 2012;Scherler et al., 2013) and other ecosystems (Okkonen and Kløve, 2011;Khoshkhoo et al., 2015).
However, there were large uncertainties in modeling soil freezing and thawing due to the complexity of phase change and coupled processes.To reduce uncertainties in modeling, uncertainty analysis method was always introduced by combining experimental data with numerical models in calibration of the models for better representativeness of reality.The generalized likelihood uncertainty estimation (GLUE) technique (Beven and Binley, 1992) is the commonly used method for uncertainty analysis in environmental modeling.Instead of searching for an optimal parameter set, the GLUE method generates ensembles of parameter sets that show equally good performance in simulations, called 'equifinality' by Beven (2006).GLUE was performed by randomly sampling the parameter space within their ranges using Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-466Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 8 October 2018 c Author(s) 2018.CC BY 4.0 License.
Monte-Carlo sampling method, and then by selecting behavioral simulations using criteria applied to performance metrics.
In this study, we performed experiments on water, heat and salt transport at two seasonal frost sites located in northern part of China.They are different in climate and in soil conditions, but are both important agricultural regions in northern part of China.The Hetao Irrigation District in China is a typical arid agricultural region suffering from soil salinization due to saline water irrigation, extensive evaporation, as well as soil freezing/thawing (Li et al., 2012).The Songyuan Irrigation District is a typical paddy rice grown region in northeastern part of China, suffering from high salinity due to overdevelopment of salinized field into agricultural field (Liu et al., 2001).Soils in both regions go through freezing/thawing during winter and suffer from salinization in spring.These two sites are crucial in water resources management of China under the concept of water-saving agriculture.Wu et al. (2016b) performed calibration on soil water and heat transport based on one plot in the experimental field in Hetao Irrigation District in Inner Mongolia and found that the influences of salt on soil freezing should be taken into account.Wang et al. (2016) conducted field experiments and analyzed water and solutes transport characteristics at these two above-mentioned sites and demonstrated that salt transport in frozen soils is more complicated than in unfrozen soils due to diffusion and solute rejection.Thus, we developed CoupModel by considering impacts of salt on freezing, and applied the new model to the agricultural sites for modeling water, heat and salt in two seasonal frost soils.The main objective was to 1) develop CoupModel by considering effects of salt on freezing point; 2) identify sensitivity of parameters; 3) analyze uncertainty in modeling soil hydrology in seasonal frost agricultural soils.

Study sites
Experiments were conducted at two agricultural sites of northern China.One site is located in Qianguo Irrigation District of Songyuan, Jilin province, China (lat: 45.24 o , lon: 124.60 o , hereafter referred as site ).This study site is typical for its soil texture classified as clay, which has a high bulk density, low porosity, and low hydraulic conductivity (Table 1).Soil profile at site Qianguo is homogeneous, with porosity of 0.46 and bulk density of 1.42 g cm -3 .The water table in this area fluctuates between 1.5 and 2.0 m.Maximum frost depth at site Qianguo is 1.2 m.Six plots (2×2 m 2 for each) were selected in a paddy field, which was cultivated with paddy rice from May to October.On 2011/10/09, 20 mm NaBr solution containing 6.5 g L -1 Br -was applied to each plot to from the initial profile for Br -.Before spraying the solution, stubbles were removed from the plots and surface was ploughed to depth of 20 cm.Annual precipitation at this site is 140 mm, and annual mean air temperature is 6.4 o C (averaged for 2012~2013).Soil profile is heterogeneous, with porosity of 0.42-0.46and bulk density of 1.44-1.53g cm - 3 (averaged over 5 plots).Water table fluctuates between 1.5 and 3 m during winter.From November 4 th to 6 th of 2012, flooding irrigation (autumn irrigation) with 250 mm water was applied to the field for leaching salt that accumulated during growing season.Soil profile mean salt content (mainly NaCl) is 0.1% g g -1 for the study site, and irrigation water electrical conductivity is 0.5 mS cm -1 .Before autumn irrigation, five plots (2×2 m 2 ) were selected for experiment at different parts of the agricultural field, and ploughed to 20 cm depth.TDR probes (Model: CS605, Campell Scientific Inc.) were installed at site Qianguo to detect liquid water content.Due to difficulty in long-term maintaining of TDR system in rural regions, only daily liquid water was manually recorded with a datalogger (TDR 100; Campbell Scientific Inc.).TDR probes were calibrated in laboratory with unfrozen soil, and the precision of calibration was maintained with R 2 of 0.97.TDR probes were then inserted horizontally into the soil pit (10 m apart from the experimental plots) from 5 cm to 100 cm depth with 10 cm interval.PT100 temperature sensors were installed at the same depth as TDR probes, and the daily temperature data were collected.100 cm with 10 cm interval were collected for determining total water content and Br -content.An electric drill (5 cm in diameter, 10 cm in length) was used for sampling frozen soil for every 10 cm depth.Total water content was determined by oven-dry method.Br -content was determined by diluting 50 g wet soil into 250 mL deionized water, and measuring the electrical potential (mV) using an electrical potential meter (MP523-06).Then the electrical potential was converted into Br -concentration by a pre-calibrated relationship between Br -concentration and electrical potential (calibration R 2 =0.99).Soil temperature and liquid water content at 4 depths (5, 15, 25, and 35 cm) from site Qianguo were used to calibrate CoupModel, while soil water storage and salt storage at various depths (0-10 cm, 0-40 cm, 0-100 cm) estimated form soil profile samplings were used to validate the model.

CoupModel_v5
Model domain covered from soil surface to 6 m depth, with unit area considered.Soil profile was discretized into 16 layers, with 10 cm thickness each layer from 0 to 40 cm, 20 cm thickness from 40 cm to 2 m, and 1 m thickness from 2 m to 6 m.Input meteorological data were hourly, and model time step was set as hourly.Numerical solution of water, heat and salt transport in soils was based on forward difference method.Model performance metrics on different output variables was calculated automatically using modules implemented into CoupModel_v5.Major model processes considered in this study were described in the following sections.

Soil water processes
CoupModel solved coupled differential equations for water and heat transfer (Jansson, 2012).Water flow in the soil matrix was described by Richards equation: (1) ( ) where is water content (m 3 m -3 ); is hydraulic conductivity (m s -1 ); is matric potential (m); is vapor diffusion coefficient (m 2 s -1 ); is vapor density (m 3 m -3 ); is the bypass flow in macro pores (m s -1 ); is depth to soil surface (positive downward) (m); and is time (s).
Vapor flow (second term inside brackets on right side of Equation ( 1)) in soil was determined by vapor gradient between two layers and diffusion coefficient, adjusted by tortuosity  !"#$ (Equation (2)).
The diffusion coefficient for free air,  ! was a function of soil temperature (Equation (A1) in Table A2), and vapor density  ! was calculated from vapor pressure (Equation (A1)), which was estimated from soil matric potential and soil temperature (Equation (A1)).
Infiltration through frozen soil was estimated separately for the low-and high-flow domains (Stähli et al., 1996).CoupModel took the preferential flow in macropores into account using a bypass routine when excess water entering the soil was routed directly to the next underlying soil layer through the high-flow domain (Jansson, 2012).Infiltration into soil was determined by the soil adsorption rate adjusted by a soil matric water adsorption coefficient  !"#$% (Equation (5)).When infiltration water was larger than soil adsorption rate, bypass flow would occur.Bypass flow in macro pores was determined by where  !"# is soil adsorption rate (m s -1 );  !"#$% is soil matric water adsorption coefficient,  ! is a geometry coefficient to describe thickness ratio to horizontal scale of each soil layer;  !"# is matric maximum hydraulic conductivity (m s -1 ); and  is pF value of soil.
Flow in low-flow domain obeyed Darcy's law and soil water retention curve was determined by Brooks and Corey (1964) equation (Equation (A3)), with the air entry at different layers to be adjusted in calibration.Hydraulic conductivity in low-flow domain was calculated from the Mualem (1976) equation (Equation (A4)).In frozen soil, hydraulic conductivity was modified for high-flow (Stähli et al., 1996).
In high-flow domain, water flow was modeled by gravitational flow under unit gradient, and hydraulic conductivity was adjusted by using impedance factor  !,! in high-flow domain (Equation (6)): impedance factor;  !"! is total water content in high-and low-flow domains (m 3 m -3 );  !" (= ! !"#$ , Equation ( 18)) and θ i ( ) are the liquid water and ice content, respectively, in the low-flow domain (m 3 m -3 ).
The hydraulic conductivity changed at the freezing front under partially frozen conditions.To prevent excessive water redistribution towards the freezing front, the hydraulic conductivity of partially frozen layers was adjusted by considering ice content influences on water flow using a factor  !" (Equation (7)).where fi c is impedance factor; and Q is heat quality, as a ratio of ice content to total water content.
Meanwhile, the influence of soil temperature on soil hydraulic conductivity was considered, using a linear increase factor  !!! and a minimum conductivity  !"!"# to adjust hydraulic conductivity at 20 o C (Equation (A7)).
Surface ponding of water may occur if the soil infiltration capacity was exceeded, otherwise the infiltration rate was equal to precipitation and rates of snowmelt.If infiltration capacity was exceeded, excess water will be transferred to the surface pool.The overland flow from surface pool was estimated by the difference between surface water storage and maximum surface pool,  !"#$ (Equation (8)).
( ) where surf a is an empirical coefficient, pool

W
is the total amount of water in the surface pool (m), and pmax w is the maximal amount of water stored on soil surface without causing surface runoff (m).
Drainage systems at two study sites were open drainage ditches, drainage at study sites was then calculated by Hooghoudt equation combined with an empirical drainage equation to constitute a manual drainage system, adjusted by initial drainage level  ! and minimum drainage level, drain spacing  !(Equation (A9)), empirical groundwater level peak value  !, and empirical groundwater flow peak value  !(Equation (A10)).Meanwhile, initial groundwater level was set for calibration.Groundwater water level was estimated by the soil saturation layer depth to surface.

Soil heat processes
Heat flow in soil was described by the heat transport equation, considering conduction, convection and latent heat flow: where is soil (containing solid, water, and ice) heat capacity (J m -3 o C -1 ); is temperature ( o C); is latent heat of freezing (J kg -1 ); is density of ice (kg m -3 ); is ice content (m 3 m -3 ); is thermal conductivity soil (W m -1 o C -1 ); is water flux (m s -1 ); is latent heat of vaporization (J kg -1 ); and is vapor flux (m s -1 ).
Upper boundary for soil heat flow was soil temperature at surface, calculated by the energy balance scheme described in Section 3.4.Lower boundary for soil heat flow was controlled by soil temperature fluctuation at 6 m depth, which was estimated using an analytical solution for soil heat conduction.
Soil thermal conductivity for both frozen and unfrozen soils was calculated from the Ballard & Arp equation (Balland and Arp, 2005), adjusted by three empirical coefficients , , and  (Equation (A11)).
Thermal conductivity from the top frozen soil layer was then corrected by using a damping function, adjusted by the maximum damping coefficient  !" (Equation (A12)).When infiltration water passed the high-flow domain, it would refreeze due to low soil temperature in frozen soils.Meanwhile, latent heat released from refreezing would melt water in high-flow domain.This would lead to redistribution of water between low-flow and high-flow domains.CoupModel considered the water redistribution and adjusted it by a heat transfer coefficient  !(Equation (A13)).

Salt tracer processes
Salt in CoupModel was simulated as a tracer migrating with water, neglecting diffusion.Salt transport was simulated as Cl -transport in soil for estimate of salt tracer flux.Salt balance in soil is calculated as: where is concentration of Cl -(kg m -3 );  !"#$% is salt deposition concentration (kg m -3 );  !"# is water flux (m s -1 ),  !"#$%% is bypass flow (m s -1 ).
Soil salt concentration for each layer was then calculated by where  !" is salt amount at each soil layer (kg m -2 );  !"# is salt adsorption rate;  is soil water content at each layer (m 3 m -3 ); Δ is soil layer thickness (m).
Salt at surface was balanced by salt in precipitation and irrigation, as well as salt loss from surface runoff.Lower and lateral boundaries for salt transport were salt leaching to groundwater, which was proportional to drainage rate.Initial salt concentration  !" , precipitation salt concentration  !"#$% , irrigation salt concentration  !"#$$#% , as well as salt adsorption coefficient  !"# at different depths were set as calibration parameters.At site Qianguo, Br -transport was converted to Cl -transport in the simulation in the validation of salt storage, Cl -storage at site Qianguo was then converted to Br -storage in comparison with field observations of Br -storage.

Energy balance processes
Surface temperature and evaporation was calculated using energy balance method, with net short-wave radiation balanced by latent heat, sensible heat and soil heat flux at surface: where v s L E is the sum of latent heat flux (J m -2 s -1 ); s H is sensible heat flux (J m -2 s -1 ) and h q is heat flux to the soil (J m -2 s -1 ).
Latent heat was calculated as below, where as r is the aerodynamic resistance (m -1 ); surf e is the vapor pressure at the soil surface (Pa or in m water); a e is the actual vapor pressure in the air (Pa or in m water); a ρ is the air density (kg m -3 ); p c is the heat capacity of air (J kg -1 o C -1 ); v L is the latent heat of vaporization (J kg -1 ) and γ is the psychometric constant.
Sensible heat was calculated as where s T is the soil surface temperature ( o C); a T is the air temperature ( o C); and as r , a ρ , p c are the same as Equation ( 13).
Soil surface heat flow was then calculated, where h k is the thermal conductivity of the topsoil layer (W m -1 o C -1 ); s T is the soil surface temperature Lq is the latent water vapor flow from soil surface to the central point of the uppermost soil layer (J m -2 s -1 ).
Surface temperature was then adjusted to make Equation (5) balanced by different fluxes at surface.Soil surface vapor pressure was determined by soil surface temperature, water potential at top layer and soil water gradient between soil surface and top layer.This was further corrected by an empirical factor, which was adjusted by an adjustment coefficient  !" , and the surface water balance (Equation (A14)), which was adjusted by maximum soil surface water deficit  !"# and maximum soil surface water excess Aerodynamic resistance for stable atmosphere was calculated using the Richardson equation.Then the aerodynamic resistance for stable atmosphere was adjusted by the momentum roughness length of soil and snow surface  !! ( !!,!"#$ ) (Equation (A16)), and the heat roughness length of surface  !! was derived from  !! and  !! (Equation (A17)).In addition, when surface was at extreme stability conditions, aerodynamic resistance was then adjusted by using a windless exchange coefficient  !,!"# !! (Equation (A18)).
Soil evaporation was adjusted by maximum soil water condensation rate  !"#,!"#$ considering the influences of condensation of water on evaporation (Equation (A19)).Net radiation was estimated by Konzelmann equation with two formulae to calculate longwave radiation and was adjusted by an empirical coefficient  !! (Equation (A20)).Snow melting was determined by solving energy balance equation in snowpack using the same scheme as soil surface energy balance calculation.Snow mass balance was then estimated based on temperature change in snowpack as well as snow age.Snow thermal conductivity was calculated from snow density with an adjustment factor  !(Equation (A21)).Soil albedo was determined by albedo of dry and wet soils, adjusted by an empirical coefficient  !(Equation (A22)).Snow albedo was determined by snow age, as well as cumulative air temperature since the latest snowfall, adjusted by the minimum snow albedo  !"# (Equation (A23)).

Soil freezing point depression function development
To solve the coupled water and heat flow equations, we needed a relation between soil temperature and soil liquid water, i.e. soil freezing characteristics.In frozen soil, when soil temperature was below zero, latent heat changed due to ice formation.When soil temperature continued decreasing, sensible heat also changed.In CoupModel, we assumed that soil is totally frozen when temperature was below  !(-5 o C), when soil temperature was between 0 and  !, sensible heat in soil was calculated as: where  is total heat stored in soil (J);  ! is latent heat of freezing (J kg -1 );  is water stored in soil (kg); ∆ is soil thickness (m);  ! is a factor accounting for the fraction of unfrozen water to soil wilting point water content;  !"#$ is the wilting point water content when the pF value of soil water is 4.2 (m 3 m -3 );  ! is density of water (kg m -3 );  ! is energy when soil is totally frozen ( ! !−  ! !"# , i.e. when soil temperature is  !,  ! is heat capacity of frozen soil, J kg -1 o C -1 );  is freezing point depression.
In modeling of soil frost, when soil was totally frozen at -5 o C, the liquid water content was determined by wilting point of soil (∆ ! !"#$  ! ), and adjusted by a coefficient  !, as depicted in Equation ( 16).Ice content in soil was calculated as: where  is total energy stored in soil (J);  is total energy stored in soil (= !, J);  ! is latent heat of freezing (J kg -1 ); ∆ is soil thickness (m);  ! is a factor accounting for the fraction of unfrozen water to soil wilting point water content;  !"#$ is the wilting point water content when the pF value of soil water is 4.2 (m 3 m -3 );  !"# is density of ice (kg m -3 ).
In CoupModel, the freezing-point depression was related to soil heat storage as below: (18) where , are empirical constants; is the pore size distribution index;  !"# is water available for freezing, kg, i.e.The second method was to relate to osmotic potential (Equation (19)).According to Banin and Anderson (1974), the relationship between freezing point and salt solution could be written as below: ( where is the freezing point ( o C); is osmotic potential (in unit cm); is a scale factor for considering the influences of salt types on the relationship (range from -2 to 2); -4 is a constant for converting osmotic potential unit from cm to MPa.
Soil salt and soil heat and water transport as well as soil freezing/thawing was connected by Equation ( 19) with osmotic potential π, and freezing point would change as soil temperature and soil salt concentration changed during simulation, since osmotic potential was determined by both soil temperature and salt concentration: where R is gas constant; T is soil temperature (K); c Cl is salt concentration (kg m -3 ); M Cl is mole mass of Cl (35.5 g mol -1 ).

Calibration approach
The sensitivity analysis and model calibration procedures are summarized in Fig. 2. We selected 58 parameters that have either shown a large influence on modeled water and heat dynamics in previous studies and/or are known to be important in sensitivity analysis (cf.Gustafsson et al., 2001;Wu et al., 2011;Metzger et al., 2015).The 58 parameters represented the major processes related to soil water, heat, radiation as well as salt transport, 19 related to soil water process, 8 related to soil heat process, 19 related to soil salt process, and 12 related to energy balance process (Table A1).We noted that 58 parameters made the calibration very inefficient, since some of the parameters were assigned to different layers and some were not so sensitive in comparison with others.We thus conducted a two-step calibration, with the first step to find out the most important parameters from different model processes based on sensitivity analysis, and the second step to calibrate the important parameters.
In the first step, the 58 parameters were tested for each site with 70000 simulations based on Monte Carlo sampling method.Each of the simulations was run with randomly selected parameter values, thus creating 70000 realizations.The most sensitive model parameters were then identified for each site based on their relative importance on performance metrics (e.g.R 2 , determination coefficient between simulation and observations, and ME, the mean deviations between simulation and observations).This was done by using the LGM (Lindeman, Gold and Merenda) method (Lindeman et al., 1980) that averages the sequential sums of squares over all orderings of regressors, which calculates the relative importance of each parameter on model performance metrics and ranks them.Based on the ranking of parameters, 8 to 11 sensitive parameters (i.e. 3 common parameters for two sites, another 5 for site Qianguo and another 8 for site Yonglian) were then selected in the second step with 10000 simulations for each site.It is important to note that the sensitive parameters may be different from site to site depending on site-specific characteristics, although initial parameters and their ranges were equivalent for all sites.From the 10000 simulations of the second calibration step, ensemble of parameter sets for each site was then selected based on statistical performance metrics (determination coefficient R 2 and mean error ME) for temperature and water at several depths (Table 2).In addition to useful information about sitespecific processes and their representations in the model, the ensemble simulation results (9 for site Qianguo, 16 for site Yonglian) simulated using accepted parameter sets were used for analysis water, energy and salt balance over the simulation period.

Freezing point depression
In Fig. 3 and Fig. 4, the sensitivity of model to freezing point depression is analyzed at 5 cm depth (the same was done for other soil depths, not shown here), based on model results from one of the behavioral simulations at each site.The influences of freezing point on soil heat are obvious (Fig. 3).
When soil freezing temperature changed from 0 o C to below zero, the relationship between soil temperature and soil heat storage changed accordingly.The model performance was improved when  The relationships between soil temperature and soil heat storage at 5 cm depth were different when various  values were assigned (Fig. 4).This indicated that different types of salt also influence soil freezing/thawing.

LGM relative importance in calibration parameters
At site Qianguo, for soil temperature R 2 at 4 depths (Fig. 5 a)-d)),  !",!"#$ (momentum roughness length of snow), which is to estimate surface aerodynamic resistance, was found to be most important.
This parameter would influence surface energy balance and eventually impact heat transport in soil profile.The other important parameter for soil temperature R 2 at 15, 25 and 35 cm depths was  !" , which is a parameter to adjust thermal conductivity of surface frozen layer.
For liquid water content R 2 at 4 depths (Fig. 5 e)-h)), the most important parameters were  !!,!"#$ , ,  ! and  ! of different depths that were related to energy balance and soil heat and water transport.
!!,!"#$ was already detected to be important for soil temperature.This indicated that surface energy balance in snowpack at site Qianguo is important for both soil heat and water transport.For soil temperature ME at site Qianguo (Fig. 5 i)-l)), important parameters were similar to soil 422 temperature R 2 , except at 5 cm depth, with  !showing the greatest importance (Fig. 5   Soil salt related parameters did not show that great importance (around 1% relative importance) to soil temperature and soil water at two sites.Even we have developed a new relationship between osmotic potential and soil freezing temperature and it has shown to be able to improve model performance on simulation of soil temperature, the parameter  did not show great importance at two sites.This indicated that for two sites,  could be assigned as fixed values for different types of salt.We only noticed the adsorption coefficients of salt at various layers show some importance to soil temperature and water.This was because they determine the osmotic potential of soil water and thus impact soil heat and water transport simulation.

Prior and posterior parameters
In Table 3, the important parameters and their posterior ranges at two sites are depicted.The posterior mean value of  !,!"# !! was reduced to 1/3 of prior mean value, and mean value for  ! was also reduced from prior at site Qianguo.The R ratio (range ratio, defined as the posterior range width ratio to prior range width) for  !,!"# !! and  ! was 0.12 and 0.18, respectively for site Yonglian (Table 3).Parameters  !!,!"#$ ,  !"!"# and  !" at site Yonglian had larger posterior mean values than prior, and the R ratio was 0.54 and 0.68, respectively.Parameters  !(1) to  !(3) at site Yonglian also obtained generally larger posterior mean values after calibration, with R ratio of 0.18 to 0.50.For  !" at site Yonglian, the R ratio was 0.22, much smaller than at site Qianguo (0.52).R ratio of  ! at 466 site Yonglian was 0.39, larger than at site Qianguo (0.18).Parameter  !,!"# !! at site Yonglian showed 467 totally different R ratio and posterior mean values from site Qianguo, with R ratio of 0.58 and posterior mean 468 value of 0.02, respectively.As discussed above,  !" ,  ! and  !,!"# !!
were important parameters at both 469 sites.They controlled soil heat and energy balance and can influence soil freezing/thawing.Salt 470 adsorption coefficient  !"# at four depths did not show large changes between posterior and prior 471 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-466Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 8 October 2018 c Author(s) 2018.CC BY 4.0 License.distributions, with R ratio from 0.74 to 0.93.The R ratio of  !from four depths at site Yonglian varied from 0.36 to 0.81.The large differences in posterior ranges of these parameters indicated these two sites have different surface water and energy balance situations.Site Qianguo is more humid in winter and has more snow events, while site Yonglian has very dry winter but more salt influences on freezing/thawing due to higher salinity at this site.At site Qianguo, parameters such as  !!,!"#$ and  !"#$% also showed to be important, and  !"# at various depths at site Yonglian were shown to be very sensitive.

Soil temperature and water
At site Qianguo, soil temperature and soil liquid water content were measured manually at daily resolution due to difficulties in installing automatic measurement instruments at farmers' land.Accepted simulations generally can capture soil temperature and water dynamics and can cover the observations within their ranges (Fig. 7 a)-b)).After calibration, the mean value of R 2 for soil temperature and soil liquid water content at 5 cm depth was 0.87 and 0.31, respectively.Meanwhile, the mean values of ME for soil temperature and soil liquid water content at 5 cm depth was -0.41 o C and -4.89%, respectively.
At site Yonglian, hourly soil temperature was obtained in calibration, and achieved high R 2 for soil temperature, with mean R 2 of 0.90 at 5 cm depth.However, soil temperature was underestimated from end of November to middle of January (Fig. 7 c)).This was mainly due to ice coverage at site Yonglian during this period.After flooding irrigation at the beginning of November, water ponding in the field (~10 cm water) was rapidly frozen and kept covering soil surface until middle of January.For this period, ice coverage disturbed water and energy balance at the site Yonglian.Even the snowpack was considered at site Yonglian and a detailed scheme for snow water and energy balance was illustrated in CoupModel, it obviously cannot describe ice coverage in our case.In the future development of the CoupModel, we recommended inclusion of a new scheme for water and energy balance on ice coverage.Due to failure of TDR in measuring water in salinized frozen soil, we only sampled total water content at site Yonglian.Soil total water content at 5 cm depth had good performance with mean R 2 from accepted simulations of 0.80.Even with only 14 sampling dates from each of five plots were selected, the calibrated model can capture soil water dynamics well.Nevertheless, we also noticed large variations in measured total water content from different plots, as indicated by error bars in Fig. 7 d).In the future development of soil water measurement methods, it is necessary to introduce more accurate measurement methods for soil water in saline frozen soils in order to obtain consistent observations of soil water dynamics during winter.

Model validation on water and salt storage
Comparison of simulated water storage with measured water storage at different soil depths is depicted in Fig. 8. Results indicated that CoupModel could predict water process well in upper 40 cm soil layer, but some large deviations mainly occurred for 40 to 100 cm at two sites between simulated and observed soil water storages.This was because the accepted simulations was derived by constraining model performance for variables (soil temperature and soil water) in upper 40 cm soil layer, and the data from 40-100 cm depth was not used for calibration.This indicated that there might be some unforeseen processes in the lower soil layer, and they would influence water processes in whole soil profile (from surface to groundwater).Since the calibration was focusing on the surface water and energy balance, and the upper layer water process was shown to be well-represented by the model, the more detailed consideration of lower layer water processes exceeded the scope of this study.Further work would include calibration of the model in the whole soil profile with more detailed measurements.The overestimation of Cl -/Br -storage at various depths indicated that the upward movement of salt with water was over-estimated.This might be due to the neglecting of diffusion and expulsion of salt in model.Cary and Mayland (1972) have shown that, the diffusion and expulsion processes in frozen soil actually played important roles in salt transport even though the convection was the major process.This was because when soil was frozen, soil solution concentrated.The concentration of soil solution would increase salt concentration gradient between soil layers.In addition, high salt concentration at low temperature would cause salt expulsion from solution due to low salt saturation (Wang et al., 2016).
However, it was difficult to measure the diffusion and expulsion of salt in frozen soil.More detailed experiments on diffusion and expulsion of salt are necessary in study of water, heat and salt coupled transport in frozen soils.Validation of soil water and salt storage more data on salt transport as well as water transport would be of importance in calibration of model, since the water and salt transport processes are tightly coupled.
where a f is the fraction of air-filled pores, 0 D is the diffusion coefficient in free air (m 2 s -1 ) and vapb d is a parameter accounting for tortuosity and the enhancement of vapor transfer observed in measurements compared with theory,  is matric potential (m, with positive value),  !"#$% is mole mass of water (18 g mol -1 ),  is gas constant,  ! is saturated vapor pressure (m),  is soil temperature (K).

Soil vapor diffusion coefficient A2
where surf a is an empirical coefficient, pool W is the total amount of water in the surface pool (m), and pmax w is the maximal amount of water stored on soil surface without causing surface runoff (m).

Surface runoff estimation A9
where  !! ,  !! are saturated hydraulic conductivities above and below water level (m s -1 ), respectively;  ! is soil depth below drainage level (m),  ! is drainage depth to soil surface (m),  ! is distance between two drainage tubes.
Hooghoudt drainage equation A10 ( ) ( ) max 0, max 0, sat sat gr z z z z q q q z z where 1 q , 2 q are maximum and minimum empirical groundwater flow (m s - 1 ), respectively; 1 z and 2 z are highest and lowest empirical groundwater level (m), respectively.

=
where s e is the vapor pressure (m) at saturation at soil surface temperature s T ( o C), ψ is the soil water tension (m) and g is the gravitational constant (g m -2 s -1 ), R is the gas constant (J o C -1 mol -1 ),  !"#$% is the molar mass of water (18 g mol -1 ) and corr e is the empirical correction factor, eg ψ is a parameter and surf δ is a calculated mass balance at the soil surface (m), which is allowed to vary between the parameters def s and excess s given in m of water.
where pool

W
is the surface water pool (m), in q is the infiltration rate (m s -1 ), s E is the evaporation rate (m s -1 ),  !"#$ is drip irrigation rate (m s -1 ), , v s q , is the vapor flow from soil surface to the central point of the uppermost soil layer (m s -1 ), def s is maximal surface water deficit (m) and excess s is maximal surface water excess (m).

Fig. 1
Fig. 1 Locations of the study sites.Site Yonlian is located at the middle part of Hetao Irrigation District of Inner Mongolia Autonomous Region, northern part of China, site Qianguo is located at Songyuan County of Jilin Province, northeastern part of China.
Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-466Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 8 October 2018 c Author(s) 2018.CC BY 4.0 License.In saline frozen soil, ice formation does not start at 0 o C, but below 0 o C. Freezing point (Equation (19)) is a parameter related to soil type, salt type and salt content.In CoupModel, was assumed as 0 o C, which was not suitable for saline soils.In this study, two methods were implemented to consider salt influences on freezing point depression.The first one was to set freezing point as a parameter in the model, and this parameter could be determined by experiments on freezing point of different saline soils.

Fig. 2
Fig. 2 Model sensitivity analysis and calibration procedures.Step 1 denotes procedures for selecting important parameters with LGM method, step 2 represents model calibration and selection of acceptable simulations.
freezing point depression was related to soil salt.Mean error (ME) for soil temperature at 5 cm depth decreased from 1.25 o C to 0.29 o C (improved by 77%) at site Qianguo, and from 2.54 o C to 1.83 o C (improved by 28%) at site Yonglian, when freezing point decreased from 0 o C to -3 o C.

Fig. 3
Fig. 3 Relationships between soil temperature and soil heat storage derived from modeling for different freezing points at 5 cm depth at a) site Qianguo and b) site Yonglian.
o C (improved by 16%) when  was 1.This indicated that in determination of freezing point, not only the salt content, but also the salt type should be taken into consideration for reducing uncertainty in modeling soil temperature.

Fig. 4
Fig. 4 Relationships between soil temperature and soil heat storage derived from modeling for different sc values at 5 cm depth at a) site Qianguo and b) site Yonglian.

Fig. 5
Fig. 5 Relative importance for parameters to R 2 a)-h) and ME i)-p) of soil temperature and soil water at 4 419 depths (5, 15, 25, 35 cm) at site Qianguo.Parameters shown in each sub-figure are the most important 420 parameter from each group, i.e. energy balance, soil heat, soil water, soil salt, and the other parameters.421

Fig. 6
Fig. 6 Relative importance for parameters to R 2 a)-i) and ME j)-r) of soil temperature and soil water at 4 430 depths (5, 15, 25, 35 cm) as well as groundwater at site Yonglian.Parameters shown in each sub-figure 431 are the most important parameter from each group, i.e. energy balance, soil heat, soil water, soil salt, then 432 the other parameters.433 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-466Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 8 October 2018 c Author(s) 2018.CC BY 4.0 License.
Fig. 7 Comparison of soil temperature and soil water at 5 cm depth at two sites.a)-b) at site Qianguo, c)-d) at site Yonglian.Red line for mean of behavioral simulations, grey lines for simulated mean values +/standard deviation (sigma), error bar in d) for standard deviation of observations from various plots.

Fig. 9
Fig.9shows the simulated salt storage in comparison with measured salt storage based on measured data at various soil depths.At site Qianguo, the Br -storage was generally overestimated by the model in comparison with measured Br -storage in whole soil profile.The simulated Br -storage showed larger uncertainty than the measured.Similarly, at site Yonglian, the simulated Cl -storage at various layers was larger than measured Cl -storage.

Fig. 8
Fig. 8 Comparison of simulated accumulated water storage with measured accumulated water storage at two sites at various soil depths.a) site Qianguo and b) site Yonglian.X error bar denotes standard deviation of observations from various plots, Y error bar denotes standard deviation of water storage from behavioral simulations.

Fig. 9
Fig. 9 Comparison of simulated salt (Br -and Cl -) storage with measured at two sites.a) site Qianguo and b) site Yonglian.X error bar denotes standard deviation of observations from various plots, Y error bar denotes standard deviation of salt storage from behavioral simulations.
!=  !!!! !" + (1 −  !" ) where  ! is soil surface frost adjustment coefficient ( o C -1 ),  !" is maximumThis is used to adjustment thermal Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-466Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 8 October 2018 c Author(s) 2018.CC BY 4.0 License.frost damping coefficient (-),  ! is surface soil temperature ( o transfer coefficient (W m -1 o C -1 ), z Δ is thickness of soil layer (m), T is soil temperature ( o C), f L is the latent heat of freezing ( u is the wind speed (m s -1 ) at the reference height, ref z (m), ib R is the bulk Richardson number, k is the von Karman constant and OM z and OH z are the surface roughness lengths for momentum and heat, respectively (mmax,cond e is maximum condensation rate (m s -1 ) for upmost soil layer to maintain water balance,  !"#$ is bare soil fraction.Soil evaporation limiting factor Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-466Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 8 October 2018 c Author(s) 2018.CC BY 4.0 License.pressure (m),  ! is air temperature ( o C) c n is cloud fraction;

Table 1
Soil physical and chemical properties at two study sites Qianguo and Yonglian.

Table 2
Criteria applied to model performance metrics in selection of behavioral simulations.
i)).  ! is for 423 estimate of snow thermal conductivity, and determines energy balance in snowpack.Site Qianguo was 424

Table 3
Important parameters and their accepted ranges at site Qianguo and Yonglian.462 a Rratio ratio of posterior parameter range to prior parameter range 463 b value without brackets is for site Qianguo, value with brackets is for site Yonglian 464 465
where dry a is albedo of dry soil, wet a is albedo of wet soil, and a k is a transform coefficient from wet to dry soil.where min a is minimum albedo of snow,  !"# is snow age (d), 1 a , 2 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-466Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 8 October 2018 c Author(s) 2018.CC BY 4.0 License.
a , and 3 a are parameter.