Simulations of water , heat , and solute transport in partially frozen soils 1

7 Experiments for soil freezing/thawing were conducted in two seasonally frozen agricultural 8 fields in northern China during 2011/2012 and 2012/2013 wintertime, respectively. Mass 9 balance was checked based on measured data at various depths. Simulation work was 10 conducted by combining CoupModel with Monte-Carlo sampling method to achieve 11 parameter sets with equally good performance. Uncertainties existed in both measurements 12 and model due to complexity in freezing/thawing processes as well as in surface energy 13 partitioning. Parameters related to surface radiation and soil frost were strongly constrained 14 with datasets available in two sites combining multi-criterion on outputs. Simulated soil heat 15 processes were better described than soil water processes given the data obtained for 16 calibration. Model performance was improved with consideration of solute effects on 17 freezing point depression. More detailed solute transport processes in CoupModel needed to 18 be improved by taking more processes such as diffusion and expulsion into consideration 19 based on more precise experimental results, to reduce uncertainty in model. Generally, 20 combination of measurement with process-based model and Monte-Carlo sampling method 21 provided an approach for understanding of solute transport as well as its influences on soil 22 freezing/thawing in cold arid agricultural regions. Incorporating more detailed descriptions of 23 processes for frozen soil in the model can be justified if uncertainties in measurements can be 24 reduced by introducing of high-precision novel technologies. 25 26


Introduction
Soil freezing and thawing processes has long been recognized for its importance in not only engineering applications (e.g., construction of roads and pipelines) (Hansson et al., 2004;Wettlaufer and Worster, 2006;Jones, 1981), but also environmental issues (e.g., soil erosion, flooding, and pollutants migration) (Seyfried and Murdock, 1997;Andersland et al., 1996;Baker and Spaans, 1997;McCauley et al., 2002).The investigation of soil freezing and thawing could result in a better understanding of water and solute distribution in soil (Baker and Osterkamp, 1989), frost heaving (Wettlaufer and Worster, 2006), waste disposal (McCauley et al., 2002), climate change in cold regions (Lopez et al., 2007).
Many experimental observations have been conducted since 1930s to study the freezing/thawing phenomenon in soil (Beskow, 1947;Edlefsen and Anderson, 1943;Spaans and Baker, 1996;Miller, 1980).For example, Beskow (1947) observed the accumulation of water towards freezing front and the influences of water, soil and solutes on freezing point depression, as well as the similarity between freezing and drying.Edlefsen and Anderson (1943) then formulated the relationship between soil temperature and freezing soil water potential by generalized Clausius-Clapeyron equation.Koopmans and Miller (1966) then tested the similarity between soil freezing curves and water retention curves.Burt and Williams (1976) measured hydraulic conductivity of freezing soil in laboratory, and this work was then put forward by others (Nakano et al., 1982;Horiguchi and Miller, 1983;Black and Hardenberg, 1991).At the same time, formulation of the soil freezing/thawing processes has been raised up.e.g., the use of Clausius-Clapeyron equation for representation of soil freezing equilibrium (Kay and Groenevelt, 1974;Groenevelt and Kay, 1974), the power relationship between liquid water content and soil temperature (Anderson and Tice, 1973), the capillary bundle model for soil freezing characteristics (Watanabe and Flury, 2008;Lebeau and Konrad, 2012), and the influences of solute on soil freezing/thawing (Bing and Ma, 2011;Azmatch et al., 2012;Wu et al., 2015).
Laboratory and field experiments on soil freezing/thawing processes have received more attention.Watanabe et al. (2013) conducted laboratory experiments to describe the influence of soil freezing on infiltration.Zhou et al. (2014) measured the water content and ice content in a freezing soil column with gamma ray attenuation and TDR method.Also, field experiments were conducted to study plot-scale or regional water, heat and solute transport during freezing/thawing from different aspects.Radke and Berry (1998) analyzed the influences of soil water, bulk density as well as microbial activities on soil water and solute transport in the field soil column experiments.Stahli et al. (2004) characterized the 4 However, there are large uncertainties in both experiments and models for soil freezing and thawing due to the complexity of phase change and coupled processes.For example, the measurement of liquid water content in frozen soil, the sampling of frozen soil, the measurements of hydraulic and thermal properties for frozen soil could be difficult due to limitations in technologies and in considering all effects on soil freezing/thawing (e.g., water, heat, solutes, soil textures as well as boundary conditions).Meanwhile, the setup of model always neglected some minor influences by taking the major one into consideration, e.g., the assumption of thermo-equilibrium of soil freezing, the neglecting of solute dispersion and expulsion in frozen soil, or even the neglecting of solute effects on freezing point depression, etc.All these would pose uncertainties to the study on soil freezing/thawing in natural conditions.To reduce uncertainties in both experiments and modeling, uncertainty analysis method is always used by combining experimental data with numerical model to calibrate the model for better representing reality.The generalized likelihood uncertainty estimation (GLUE) technique (Beven and Binley, 1992) is the commonly used method for uncertainty analysis in environmental modeling.
In models for soil water, heat and solute transport (e.g., CoupModel, HYDRUS, SHAW), there are many parameters related to different coupled transport processes, also including non-linear responses , especially when considering soil freezing/thawing and solute transport.The parameters in these models are possible to measure with independent methods but those are big challenges because of high variability in the environments and many temporal and spatial scale related dependencies.Also, model structures are always simplified for description of some processes.Thus, for investigation of coupled processes in seasonally frozen agricultural field, it is important to use the GLUE method combining process-based model to unveil the freezing/thawing phenomenon.
In this study, we conducted field experiments on water, heat and solute transport in two sites in northern part of China.They are special both for the climates and the soils in cold Hydrol.Earth Syst.Sci. Discuss., doi:10.5194/hess-2016-507, 2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: November 2016 c Author(s) 2016.CC-BY 3.0 License.regions in China.The coupled transport of water, heat and solute during wintertime is common for these two sites.Also, due to limited tools in measurements, only the common variables were observed using common methods during experiments.
The main objective was to search for constrain with the current available data and models by considering also explicit salinity impacts on freezing in the model.With the collected data from field, the model was combined with Monte-Carlo sampling method to 1) investigate how well simulated water and solute dynamics corresponded to measured data; 2) identify variability in modeling of water, heat and solute using ensemble simulations; and 3) discuss the influences of solute on soil freezing/thawing based modified freezing point depression function.

Study sites
Studies were conducted at two experimental fields in north China, during 2011/2012, and 2012/2013 wintertime, respectively.One site was located in Qianguo Irrigation District, Songyuan, Jilin (Site NE) (Fig. 1).Annual precipitation in Site I was 451 mm and mean monthly air temperature was 5.1 o C. The study site is typical for its soil texture of clay, which has a high bulk density, low porosity, and low hydraulic conductivity (Table 1).The water table in this area fluctuated between 1.5 and 2.0 m.Maximum frost depth in Site I was 1.2 m.
In Site I, six plots (2×2 m 2 , denoted as P1 to P6) were selected in an agricultural field, which was cultivated with rice from May to October.On 2011/10/09, 20 mm NaBr solution containing 6.5 g L -1 Br -was applied to six plots to form the initial profile for Br -.Before spraying of the solution, stubbles were removed from the plots and surface was ploughed to depth of 20 cm.Field experiment in Site I was conducted during 2011/2012 wintertime.
The other site was located in Hetao Irrigation District, Inner Mongolia, China (Site IM) (Fig. 1).Annual precipitation in this site was 140 mm.Annual mean air temperature 6.4 o C. Soil texture in this site is characterized as silt loam, with porosity of 0.42~0.46and saturated hydraulic conductivity of 3.84×10 -5 m s -1 .Water table was kept between 1.5 and 3 m for the winter time, and three irrigation events occur every year in May, July, and November.Soil salt content (mainly NaCl) was 0.1% g g -1 for the study field, and irrigation water with electrical conductivity of 0.5 mS cm -1 .Before autumn irrigation, five plots (2×2 m 2 , denoted as D1-D5) were selected at different parts of the agricultural field, and ploughed to 20 cm depth.Field experiment in this site was conducted during 2012/2013 wintertime.

Experimental design
TDR probes with 15-cm long and three-rod (CS605) were installed in Site NE to detect liquid water content.A datalogger (TDR 100; Campbell Scientific Inc.) was connected to the probes and recorded daily water data.TDR probes were inserted horizontally into the soil pit (10 m apart from the plots) from 5 cm to 100 cm with 10 cm interval.TDR probes were calibrated in laboratory with unfrozen soil, and the precision was maintained within R 2 of 0.97.PT100 temperature sensors were installed at the same depth as TDR probes, and the daily temperature data was also collected.During soil freezing/thawing, sampling was conducted at 7 dates (2011/10/09, 2011/11/09, 2011/11/25, 2011/12/20, 2012/02/15, 2012/04/10, 2012/04/20), and soil samples from 5 to 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 soil from depth to 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 (R 2 =0.99).Total water content and Cl -content in Site IM were sampled at 14 dates (with around 25-d interval) from October 2012 to April 2013April (2012April /10/16, 2012April /10/27, 2012April /11/10, 2012April /12/04, 2012April /12/15, 2012April /12/26, 2013April /01/05, 2013April /01/14, 2013April /01/25, 2013April /03/05, 2013April /03/14, 2013April /03/25, 2013April /04/07, 2013/04/18)/04/18).The sampling and measurement methods for total water content and Cl -content were the same with those in Site NE.Hourly soil temperatures at 5, 15, 25 and 35 cm depth were recorded by the PT100 temperature sensors from the micro-meteorological station in the field.Groundwater table depth was measured for every 5 d, and the frequency was increased to every 1 d during the autumn irrigation period (2012/11/4 to 2012/11/15).Meteorological data for two sites, e.g., air temperature, humidity, radiation, wind speed, and precipitation, were obtained from the nearest meteorological station with hourly-resolution.

CoupModel theory
In seasonally frozen soil, the transport of water, heat, and solute in soil profile is coupled with lower boundary (groundwater) and upper boundary (atmosphere) (Fig. 2).Water transport processes in CoupModel could be described by combining Darcy's law with mass conservation law: Hydrol.Earth Syst.Sci. Discuss., doi:10.5194/hess-2016-507, 2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: November 2016 c Author(s) 2016.CC-BY 3.0 License.

 
where  is water content (m 3 m -3 ); is the bypass flow in macro pores (m s -1 ); z is depth to soil surface (positive downward) (m); and t is time (s).
Hydraulic conductivity in freezing/thawing soil is modified by dividing water flow domain into high flow and low flow domains (Stähli et al., 1996), and adjusted by using impedance factors, respectively.
Heat flow in soil is described by the heat transport equation, considering conduction, convection and latent heat flow: where C is soil (containing solid, water, and ice) heat capacity (J m -3 o C -1 ); T is temperature w q is water flux (m s -1 ); v L is latent heat of vaporization (J kg -1 ); and v q is vapor flux (m s -1 ).
Thermal conductivity for frozen/unfrozen soil is calculated by the Balland and Arp (2005) method considering the influences of soil components on heat flow.
Solute in CoupModel is considered to transport with water, neglecting diffusion.Solute transport is converted into Cl -transport in soil: where Cl c is concentration of Cl -(kg m -3 ); and w q is water flux (m s -1 ).
Upper boundary for model is atmosphere and snow layer is also taken into consideration.Lower boundary is saturated soil layer and drainage is calculated by the combination of empirical drainage equation with Hooghoudt drainage equation The detailed descriptions of water, heat, and solute transport processes as well as model boundaries could be found in the CoupModel manual (Jansson and Karlberg, 2004).Equations used for this simulation work were detailed listed in Table S2 in Appendix.

Modification of freezing point depression functions
In CoupModel, the freezing-point depression (Fig. 3(a)) is described as below: Hydrol.Earth Syst.Sci. Discuss., doi:10.5194/hess-2016-507, 2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: November 2016 c Author(s) 2016.CC-BY 3.0 License.T is introduced with values below 0 (from -3 to 0 o C), the freezing point depression will be like in Fig. 3(b).
When calibrating the model in terms of freezing/thawing, parameter 0 T needs to be determined with respect to salt influence on soil freezing/thawing.
As the influence of salt on freezing point is mainly dependent on the osmotic potential of soil solution, a third relationship between freezing point and osmotic potential is built.
According to Banin and Anderson (1974), the relationship between freezing point and salt solution could be written as below: where c is salt concentration (mol L -1 ); m N is number of ions to which the salt molecule dissociates; Z is valency of salt.
In CoupModel, osmotic potential ( (cm)) is a function of salt concentration where R is gas constant (8.31 J mol -1 K -1 ); T is temperature ( o C); c is salt concentration (mol L -1 ).
Substituting Equation (6) to Equation (5) will obtain: For simplicity, Equation (7) could be expressed as below with a linear relationship between freezing point and osmotic potential:

Uncertainty analysis approach 1) Model performance indicator selection
The likelihood is necessary for a formal objective selection of accepted simulations based on the uncertainty of also the measurement method.However, this requires information about various errors in the measurement approach, which is tricky to obtain without making very extensive investigations.There are many different kinds of substitutes for likelihood functions (Beven and Freer, 2001) when evaluating model performance, and they are widely used in hydrology (e.g., Li et al., 2010;Besalatpour et al., 2012;Pathak et al., 2012).Here, the Nash-Sutcliff index was selected as a performance indicator to be used to reject nonbehavioral simulations.This function is calculated as follows:  is the indicator for each model run with parameter set  , N is the total number of measurements, i Y is the measured value for the t h i measurement and ˆi Y is the corresponding output of the model, Y is the average of the measured data.If the model predicts the measurements perfectly, we have î has the same goodness of fit as using the average of the measured data for every situation.

2) Model calibration with Monte-Carlo sampling method
In this work, we calibrated the model in two study sites using data from one winter period within the GLUE framework, respectively.Since this study is to discuss the model performance and parameter uncertainty for a process-based model in simulation of water, heat and solute transport in frozen soils, we would focus more on the calibration of the model instead of validation.
In a physically based model assuming possible variability between soil layers, around 100 parameters are used for simulation.However, for most parameters, they are not necessarily to be calibrated according to the interest of the modeling work.Thus, in this study, a set of parameters that show high sensitivities were selected for calibration, as shown in Table S1.These parameters were selected on the basis of one-parameter-at-a-time sensitivity analysis, which was conducted before setting up the model.For the other parameters, they were set as fixed values based on experimental results or references to previous research (Wu et al., 2011a;Gustafsson et al., 2001;Metzger et al., 2015).Then, a uniform prior distribution was assigned to each parameter and 70,000 random parameter sets The simulated variables for two sites (temperature, liquid water content at 5, 15, 25, and 35 cm depth for Site NE; temperature, total water content at 5, 15, 25, and 35 cm depth, and groundwater level for Site IM) were compared with mean of measured from multi-plot, and NSE R 2 was used to constrain model performance to achieve the accepted simulations.At Site NE, NSE R 2 >0.7 for soil temperature at multi-depth and NSE R 2 >-1 for liquid water content at multi-depth were chosen to constrain the simulations.At Site IM, NSE R 2 >0.8, NSE R 2 >-0.5, and NSE R 2 >-0.5 were used for constraining temperature, total water content and groundwater table depth, respectively.These criteria were chosen based on the cumulative distribution of the model performance, by reducing the number of accepted simulations to around 200, and keeping the model performance (NSE R 2 ) of accepted simulations higher than 80% of all the simulations.Finally, 204 and 222 accepted simulations were obtained for Site NE and Site IM, respectively.

Water and solute balance analysis
Changes of water and solute storage for different soil horizons showed typical patterns (Fig. 4).Seven sampling dates at Site NE divided the whole experimental duration into six periods.While at Site IM, 13 sub-periods were obtained by 14 sampling dates during the experiment.At Site NE, water storage tended to increase from Period 1 to 5 (Fig. 4(a)), for 0-10, 0-40, and 0-100 soil zones, except for Periods 2 and 3, when water storage tended to decrease for some zones.For Period 2, water storage decrease in 0-10 cm zone was mainly due to evaporation loss because solute storage (Fig. 4(b)) in 0-10 cm zone of this period increased.Besides, solute storage change in three zones during Period 2 were similar, which meant that water and solute change in soil layer lower than 10 cm depth did not influence water and solute storage change in 0-10 cm zone.
Similarly, in Period 3, water loss in 0-10 cm soil layer was also due to evaporation.In Periods 4 and 5, water storage increased largely in soil profile, because of upward movement of water under temperature and potential gradients during soil freezing.In Period 6, water storage in 0-100 cm depth generally decreased, because in this period, soil was thawing, evaporation, runoff would cause large amounts of water loss from soil profile.Solute storage in Periods 5 and 6 showed less change or slight decrease due to loss of water from soil profile.
Water storage at Site NE (Fig. 4(c for exception.This was because Period 2 is the period when soil was going through extensive evaporation for the fact that water storage at 0-10 cm depth decreased while solute storage (Fig. 4(d)) at the same depth increased.Then during the whole winter, water storage in 0-100 cm soil depth kept increasing.Solute storage showed continuous increase in the whole winter-time.Comparing water and storage changes in all three zones for each period, it could be found that changes occurred in the whole soil profile.It was not like that at Site NE, where changes occurred only in the upper layer for some time.This might be due to the climatic differences in two sites during experiments.At Site NE, more precipitation occurred during the winter-time of 2011/2012, which could disturb water and heat as well as solute transport in soil profile due to infiltration of water, snow cover, and re-freezing of soil.The 0-10 cm soil layer at Site NE would be influenced largely by this disturbance and thus show more frequent changes in water and solute storages during experiment.Site IM was characterized as dry winter with sparse precipitation during the experiment.Soil water and solute transport in this site was mainly influenced by the water potential and temperature gradients in soil profile.High groundwater level in this region would also enhance the upward movements of water and solute during winter-time 2012/2013.

CoupModel performance and parameter uncertainty analysis
Most parameters show similar posterior ranges as the prior ranges (Table S1).Large differences between prior and posterior ranges were detected for C md at Site NE, and for ψ eg , s def , r a,max -1 , ψ a (1), ψ a (3), ψ a (4) at Site IM.The correlations between parameters were normally low, showing that most of the parameters can be taken as independent (Table S1).
For the parameter showing large differences in prior and posterior range ratios, NP was always larger than 1, indicating there were some parameters showing high correlations with it.It should be noticed that the posterior range was obtained based on the substantial improvement of model performance even if the change of most parameters distributions were small.The correlation between parameters should be taken into consideration also when the correlation is low since the combination of many parameters still improved model performance and make interpretation of individual parameters uncertain.However, it does not mean that there must be direct relationships between two parameters with high correlation coefficients.Sometimes the "mass correlation" might exist in statistical analysis work, which means that the correlation between two parameters does not result in the conclusion that the two parameters are connected to each in the model, it is just a statistical result based on large amounts of samples (here the samples are accepted parameter sets).Thus, a further step of analysis work was conducted (not shown in the study) to plot the scattered relationships between the parameters with correlation coefficients (absolute values) larger than 0.2.Then, the plots were checked one by one based on the scattered relationships and the processes related to these parameters.Finally, 6 parameters in each site were determined as sensitive parameters in calibration, and the cumulative distributions of these parameters are illustrated in Fig. 5 and Fig. 6.For Site NE, the model-sensitive parameters were r a,max -1 , d 1 , z 0M,snow , C md , s k , r k1 , and they were connected radiation and evaporation process (r a,max -1 , z 0M,snow , s k , r k1 ) and frozen soil heat process (C md , d 1 ).For Site IM, the six sensitive parameters belonged to radiation and evaporation process (ψ eg , s def , r a,max -1 , r k1 ) and soil water process (ψ a (1), ψ a (4)), respectively.These parameters showed very different posterior distributions in comparison with the prior uniform distributions, and they showed rapid increase in cumulative probability for certain part of the pre-set ranges.
When looking into the processes that were related to these parameters, it could be found that radiation and evaporation process in both sites seemed to be more sensitive in parameterization.For example, r a,max -1 and r k1 are two parameters accounting for windless conditions and radiation estimations, respectively.The distributions of these parameters indicated that the proper choice of parameter ranges and distribution related to soil evaporation and radiation process would be of importance in obtaining high model performance.Besides, due to the site-specified conditions in simulation, the specific processes in each site should also be taken into consideration.For example, at Site NE, snow and frost processes showed more sensitive feedbacks with model performance.At Site IM, the water processes in soil and the surface energy balance needed to be taken into account more carefully in calibration of model.
In Table 2, the accepted range ratio for model performance for multi-variable is depicted, as well as the number of parameters showing strong correlation with each variable.
At Site NE, the accepted ratios for NSE R 2 with respect to all the variables of interest were well constrained to the 50% best, using the selected criteria.The model performance (NSE R 2 ) for soil temperature at four selected layers was better constrained than total water content and liquid water content at respective layer.This might firstly be due to the fact that soil temperature data was more sufficient (with daily resolution) in comparison with total water content (7 records for each layer) and liquid water content data (around 30 records for each layer).The measurement of soil temperature was automatic with higher precision in comparison with the manual measurement of soil total water content.
Besides, the measurement of liquid water content with TDR was not totally automatic in uncertainties in the data obtained by TDR (Topp et al., 1980;Stähli and Stadler, 1997).The second reason for the poorer performance in modeling water than in temperature was that the water process in frozen soil is much more complex than heat transport, with ice formation, refreezing and infiltration and drainage considered.These processes are actually influenced by each and result in a high nonlinearity in the model.
The modeling of water and temperature in seasonally frozen soil was reported to have a trade-off, because of the coupled relationships between water and heat transport processes (Wu et al., 2011).At Site NE, the higher precision in temperature data prevailed in obtaining better model performance.It would as a result cause the decrease in model performance for water when constraining the model with selected criteria and aimed size for accepted simulations.For the number of parameters showing strong correlation with model performance of each variable, it could be noted that accepted model performance for each variable was correlated to certain number of parameters in their accepted ranges.This relationship indicated that, when aiming at achieving certain level of model performance for some variables of interest, the parameters related to the processes would be constrained accordingly.They should also be taken into consideration to try to get rid of the condition that one parameter showed opposite feedbacks to variables of different processes or different layers.This would also cause the trade-off between model performance indicators in variables, especially when a multi-objective constraining study was required, both water and temperature were constrained for accepted simulations.
Similarly, at Site IM, the model performance was also well constrained to the 30% best of all the simulations.Soil temperature and total water content model performances were constrained to the 5% best, which were higher than at Site IM.This was mainly because that at Site IM, the soil temperature used for calibration was with 1-h interval recorded data, this would result in a better model performance.The application of criteria for soil temperature model performance would not reduce the number of simulations largely.Even though the total water content data was obtained manually at several dates of the experiment, with relatively large uncertainties, the constraining of total water content model performance would produce reasonable number of accepted simulations.For groundwater table depth, the 23% best of the simulations was obtained, showing good constraining of this variable.Also, several parameters were shown to be tightly related to model performance of each variable,  when soil was at the beginning of freezing or end of thawing.Most of the data was obtained when soil was deeply frozen with less water movement in upper layer.This indicated that model parameters might be well constrained for frozen conditions with less water movement.
These parameters would not be well applicable for conditions when soil was not frozen deeply with larger water movement in profile or in certain layers.It also implied that it might be good to consider the seasonal characteristics for some parameters for achieving modeling results with consistent performance for the whole period.
In this study, data for non-frozen season was not accessible, further work will be necessary to consider the seasonal patterns in calibration work.When looking into the mean of measured total water content in comparison with modeled mean, total water content at 5 cm at Site IM was generally underestimated by the model, especially for the frozen season (after December).And the range for the measured total water content at 5 cm depth reached as high as 60% cm 3 cm -3 , which indicating an over-saturation condition.Eye-inspection during sapling also verified this phenomenon that soil samples at surface layer was detected with some ice-lenses existing, as shown in Fig. S1.These ice-lenses would result in large water content in soil and also cause the frost heave in soil profile.This ice-lens at Site IM was mainly due to the intensive irrigation before soil freezing, which maintained high soil water content profile during freezing.When soil began freezing, large amount of water would move upward to occupy the soil pores, and surplus water would accumulate at freezing front, forming the ice lens (Konrad and Morgenstern, 1980;Konrad and Morgenstern, 1981;Watanabe and Mizoguchi, 2000).
Liquid water content at Site NE (Fig. 7(c)) was underestimated by the model, in comparison with the measured.This might be due to the fact that, the influence of solutes on freezing/thawing was not taken into account in the model, which would result in a lower modeled liquid water content during soil freezing/thawing.Studies by many researchers (Banin and Anderson, 1974;Stähli and Stadler, 1997;Wu et al., 2015;Wu et al., 2016) have shown that the influences of solutes on soil freezing and thawing should not be neglected in modeling water, heat and solute transport in frozen soil, because the freezing point depression effect was strong.Groundwater table depth at Site IM was also underestimated by the model (Fig. 7(f)).In the model, drainage was described by combination of empirical with physical drainage equations.The returning flow was not taken into consideration, as well as the freezing of water in drainage system during the winter.Due to the rapid freezing after irrigation, most of the drained water in the ditches were frozen immediately, and was melted until the end of soil thawing.This would result in a high water level in drainage system, which would make it impossible for efficient drainage during soil freezing/thawing periods.A combination of more detailed agricultural field drainage processes with water body freezing/thawing in cold regions would be of high interest in the future to better predict dynamics of groundwater.

Uncertainties in water-solute dynamics simulations
Comparison of simulated water storage with measured water storage in different soil profiles is depicted in Fig. 8. Results indicated that CoupModel could predict water process well in upper 40 cm soil layer, and some errors occurred for prediction of layer between 40 and 100 cm.This was because the accepted simulations was derived by constraining model performance for variables in upper 40 cm soil layer, and the data from layers lower than 40 cm was not used for validation.This indicated that there might be some other processes in lower layers there were not included in the upper 40 cm layer, and would influence water processes in whole soil profile (from surface to groundwater).Since the model calibration work 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 be conducted to calibrate the model for the whole soil profile with more detailed measurements.
Fig. 9 shows the simulated salt storage in comparison with measured salt storage based on measured data at various soil layers.For Site NE, the Br -storage was generally overestimated by the model in comparison with measured Br -storage in whole soil profile.Also, the simulated Br -storage shows larger uncertainty than the measured.Similarly, at Site IM, the simulated Cl -storage at various layers was larger than measured Cl -storage.Also, simulated Cl -storage was shown with larger uncertainty.In the calibration work, only outputs for soil heat (e.g., soil temperature) and soil water (e.g., soil water, liquid water, groundwater) processes were constrained with certain criteria.Due to the assumption of only convection for solute transport in soil, the simulated solute dynamics was tightly related to water dynamics.As is known that, during soil freezing period the upward of water and solute was the major process.The over-estimation of Br -storage at various depths indicated that the upward movement of Br -with water was over-estimated.This might be due to the neglecting of diffusion and repulsion of solute in model.Studies by Cary and Mayland (1972) have shown that, the diffusion and expulsion processes in frozen soil actually played important roles in solute transport even though the convection was the major process.This was because when soil is frozen, soil solution would condense.This would increase solute concentration gradient between frozen layer and unfrozen layer.Also, high solute concentration at low temperature would cause solute expulsion from solution due to low solute saturation.
However, it is difficult to obtain the diffusion and expulsion of solute in frozen soil.More detailed experiments on diffusion and expulsion of solute would be of high interest in study of water, heat and solute coupled transport in frozen soils.Large uncertainties in simulated solute storage in both sites indicated that the simulation of solute transport needed to be taken into account more carefully.e.g., more data on solute transport as well as water transport would be of importance in calibration of model, since the water and solute transport processes are tightly coupled in CoupModel.

Influences of solutes on soil freezing/thawing
In modeling of transport processes in saline frozen soil, the influence of solute on freezing point is also another uncertainty.In CoupModel, the formulation of freezing point depression only takes soil types into consideration, while the solute has shown to be a more important factor in freezing point depression (Banin andAnderson, 1974, Wu et al., 2015).
CoupModel takes the freezing point of soil as 0 o C, which is not reasonable in most agricultural soils because minerals commonly exist in agricultural fields.
In Fig. 10 where 1 q , 2 q , 1 z and 2 z are parameters related to hydro-geological conditions.

Empirical
where E is the total energy content in the soil, f L is the latent heat of  , ( 1) ( ) max , min () where pool W is the surface water pool, in q is the infiltration rate, s E is the evaporation rate, ,

d
are empirical constants,  is the pore size distribution index.In saline frozen soil, ice formation does not start at 0 o C, but at freezing point of 0 related to soil type, solute type and solute content.In this approach, 0 T needs to be defined as a parameter in CoupModel, and this parameter could be determined based on experiments or by calibration with soil temperature data for frozen soil.In the present version of CoupModel, 0 T is actually taken as 0 o C, as shown in Fig.3(a).In this development of the model, 0

T
is the freezing point ( o C);  is osmotic potential (cm); sc is a scale factor for considering the influences of solute types on the relationship (range from -2 to 2); -4 is a constant for converting osmotic potential unit from cm to MPa.Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-507,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: November 2016 c Author(s) 2016.CC-BY 3.0 License.
Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-507,2016   Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: November 2016 c Author(s) 2016.CC-BY 3.0 License.this study.The measurement was taken at some irregular days by connecting the TDR cables to the datalogger.The measurement of liquid water content with TDR was influenced by a lot of factors, such as soil texture, soil solute and water conditions, these would result in large Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-507,2016   Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: November 2016 c Author(s) 2016.CC-BY 3.0 License.indicating the complexity and nonlinearity in model, and the uncertainty in obtaining a multiobjective model constraining result.

Fig. 7
Fig.7shows the modeling results with selected variables in two sites.Uncertainties in measurements and modeling results were added for better understanding of uncertainties in calibration work.Soil temperature at 5 cm depth (Fig.7(a), (d)) show small uncertainties both in measurements and modeling.Total water content at 5 cm depth show relatively large uncertainties in two sites, and larger uncertainties for modeled total water content existed Hydrol.EarthSyst.Sci.Discuss., doi:10.5194/hess-2016-507,2016   Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: November 2016 c Author(s) 2016.CC-BY 3.0 License.
Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-507,2016   Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: November 2016 c Author(s) 2016.CC-BY 3.0 License.However, as could be detected from Fig.8(a), water storage was generally under-estimated for different soil layers.This was in conflict with results from Br -dynamics.Br -was added to soil surface before freezing/thawing, and the original Br -in soil could be neglected.The only changes in Br-storage in each layer could be attributed to transport to or from adjacent soil layers.At Site IM, the Cl -was used as tracer.Spatial variability in Cl -caused large uncertainty in measured Cl -storage at various depths (Fig.8(b)).
and Fig. 11, the sensitivity of model results to freezing point depression was analyzed for 5 cm soil layer, based on the proposed freezing point depression functions.The influences of freezing point on soil heat are obvious in Fig. 10.Also the model performance was improved by considering freezing point depression due to solute.Mean error (ME) for Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-507,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: November 2016 c Author(s) 2016.CC-BY 3.0 License.soil temperature 5 cm depth decreased from 1.25 o C to 0.29 o C at Site NE, and from 2.54 o C to 1.83 o C when freezing point decreased from 0 o C to -3 o C.This indicated that the neglecting of solute effects on freezing/thawing of soil would cause uncertainties in simulation results.As shown in Fig.11, the relationships between soil temperature and soil heat storage at surface 5 cm layer are different when different values of sc are assigned.ME decreased from 1.35 o C to 0.64 o C when sc changed from 0 to 1 for Site NE.For Site IM, ME decreased from 2.54 o C when sc is 0 to 2.14 o C when sc is 1.This indicated that in calculation of freezing point, not only the solute content, but also the solute type should also be taken into consideration for reducing uncertainty in modeling soil temperature when soil solute dynamics in frozen soil is not negligible.In Fig.12, model performance changes in simulating soil temperature at multi-depth were compared before and after adding solute influences on freezing point depression, by using Equation (8).sc was calibrated with ranges from -2 to 2, based on the 204 and 222 calibrated parameter sets with other parameters.The cumulative distribution for NSE R 2 is plotted in Fig.12for comparison.Model performance was improved for different soil depths in both sites, especially for Site NE.This was because at Site NE, only daily mean data was used for calibration of soil temperature.This could cause large uncertainty in simulation of frozen soil.The consideration of solute effects on soil freezing/thawing could improve model performance in simulating heat as well as water processes.This in-turn compensated the uncertainty in inputs.For Site IM, model performance in simulating soil temperature was significantly improved in upper layers (e.g., 5 cm and 15 cm depths).This indicated that detailed description in solute influences on soil freezing/thawing could reduce uncertainty in modeling, and thus improve model performance in modeling surface water and energy balance.and solute during freezing/thawing were simulated and compared with measurements in two seasonally frozen soils in northern China.The uncertainties in both measurements and model were discussed using Monte-Carlo sampling method.Seasonal patterns of water and solute dynamics were detected in both sites during soil freezing/thawing.Water in different soil layers was influenced by both freezing/thawing and evaporation.Solutes tended to accumulate in upper soil layer during the winter.Some model parameters related to radiation/evaporation, soil freezing and snow pack as well hydraulic Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-507,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: November 2016 c Author(s) 2016.CC-BY 3.0 License.properties were possible to constrained to a more narrow range based on the applied multicriterion constraining soil temperature, liquid or total water content and groundwater table.
However, trade-off existed between different criteria as well as between different model outputs.Solute has shown to be a key factor influencing soil freezing/thawing.Diffusion and expulsion of solute in frozen soils could result in uncertainties in estimating both water and solute transport with numerical model.Modification of freezing point depression function has shown to be necessary in simulation of mineral-contained soils.Uncertainty in inputs could also be compensated by more reasonable representation of freezing/thawing processes in agricultural soils.CoupModel combined with ensemble simulations method provides an approach for understanding water, heat and solute dynamics as well as parameter uncertainty in seasonally frozen soils.More phenomena like the diffusion and expulsion of solute could not be clarified by the limited amount of data.Detailed experiments on solute transport mechanism would be of high interest in investigation of salinization in cold arid agricultural districts.Future work using more precise measurements and maybe also incorporating more processes into the model are still of high interest in understanding coupled processes in seasonally frozen soil.

Fig. 3 Fig. 4
Fig. 3 Modification of the energy-temperature relationship in saline frozen soil.

Fig. 5 Fig. 6
Fig. 5 Cumulative distribution for selected sensitive parameters at Site NE.Grey dashed line denotes prior distribution, and black solid line denotes posterior distribution.

Fig. 7 Fig. 8
Fig. 7 Comparison of simulated results with measured data in two sites.(a), (b) and (c) are soil temperature, total water content and liquid water content at 5 cm depth at Site NE; (e), (f) and (g) are soil temperature, total water content at 5 cm depth and groundwater water table depth (positive upward) at Site IM.Band shows the accepted range for each variable with minimum and maximum values, and error bar shows the standard error range for all plots.

Fig. 9
Fig. 9 Comparison of simulated solute storage with measured in two sites.(a) for Site NE and (b) for Site IM.X and Y error bars denote standard errors for measured and simulated results, respectively.

Fig. 10
Fig. 10 Relationships between soil temperature and soil heat storage for different freezing point at 5 cm depth at (a) Site NE and (b) Site IM.

Fig. 12
Fig. 12 Comparison of cumulative distribution of NSE R 2 for soil temperature at different soil depths without or with consideration of solute effects on freezing point depression.(a) 5 cm, (b) 15 cm, (c) 25 cm, and (d) 35 cm depths.Red dashed and solid lines denote cumulative distribution of NSE R 2 without consideration of solute effects on freezing point depression with ensemble simulations.Blue dashed and solids lines denote cumulative distribution of NSE R 2 with consideration of solute effects on freezing point depression with ensemble simulations.
vapor pressure at the soil surface, a e is the actual vapor pressure in the air, a  is the air density, p c is the heat capacity of air, v L is the latent heat of vaporization and  is the psychometric constant.Latent heat flux (E.24)   =     (  −   )   where s T is the soil surface temperature, a T is the air temperature and Sensible heat flux Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-507,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 28 November 2016 c Author(s) 2016.CC-BY 3.0 License.middle of uppermost soil compartment temperature, 1 z  is the depth of the uppermost soil compartment and , v sL q is the latent water vapor flow from soil surface to the central point of the uppermost soil layer.vapor pressure at saturation at soil surface temperature s T ,  is the soil water tension and g is the gravitational constant, R is the gas constant, M is the molar mass of water and corr e  is a parameter and surf  is a calculated mass balance at the soil surface, which is allowed to vary between the parameters def where u is the wind speed at the reference height, ref z , i b R is the bulk Aerodynamic resistance Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-507,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 28 November 2016 c Author(s) 2016.CC-BY 3.0 License.40 Richardson number, k is the von Karman constant and OM z and OH z are the surface roughness lengths for momentum and heat, respectively.Syst.Sci.Discuss., doi:10.5194/hess-2016-507,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 28 November 2016 c Author(s) 2016.CC-BY 3.0 License.Syst.Sci.Discuss., doi:10.5194/hess-2016-507,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 28 November 2016 c Author(s) 2016.CC-BY 3.0 License.

Table 1
Soil properties for two sites.
o k is the conductivity of the organic material at the surface, s T is the surface temperature, 1 T is the temperature in the upper most soil Upper boundary condition Hydrol.inq is the water infiltration rate, s  is the dry bulk soil density.
ice w is the mass of water available for freezing and r is the freezing- 3 d are parameters and  is the pore distribution index.