Simulating Water and Pollution Exports from Soil to Stream during the Thawing Period at the Small River Basin Scale

: A physical model was developed to describe the soil-to-stream export processes of water and pollutants in a small river basin during the soil thawing period. The hydrological and pollution transport and transformation behaviors in paddy- and corn-dominated catchments were numerically simulated; the impacts of the pollution concentrations, interactions between the migrated water and pollutants in the soil, and pollutant transformations during the pollution export were coupled. Experimental ﬁeld data from the Heidingzi river basin during the soil thawing period were used to calibrate the model parameters and evaluate the performance. The mass of the dissolved pollutants from soil particles in the migrated soil pore water was the key factor affecting the pollution export into the streams; the water content directly affected the pollution export. The concentration of the pollutants peaked when the initial exported water was high. The pollutant transport processes inﬂuenced the pollution export more signiﬁcantly after the soil water was signiﬁcantly reduced. The N-S efﬁciency coefﬁcients between the simulated and monitored ﬂow rates and the pollution concentrations at the outlets of the paddy- and corn-dominated catchments were >0.60 and >0.54, respectively. The system deviations between the simulated and monitored ﬂow rates and the pollution concentrations were <10% and <15%, respectively. The proposed model effectively described the water ﬂow, pollution transport and transformation processes.


Introduction
In seasonal freeze-thaw regions, soil freezing affects the biological processes of microorganisms and the transformation of nitrogen and phosphorus. With the reduction of vegetation and roots, a lot of pollutants accumulate in the soil [1]. Because the soil and snow melt together in spring, the non-point source pollution in the cold regions is cumulative and sudden [2]. However, due to the complexity of the pollutant migration and transformation in the frozen soil thawing process [3], it is difficult to quantitatively analyze the non-point source pollution during this period. Therefore, it is of great significance to simulate the water and pollution export from soil to stream during the thawing period. In unfrozen soil, the transport process of dissolved pollutants can be described using a classic convection-diffusion equation with sink or source terms [4]. However, freezing and thawing affect the migration and transformation of pollutants, primarily in the following three aspects: (1) Freeze-thawing changes the structure of the soil aggregates, which affects the permeability of the soil [5,6]; the contact areas between the soil particles and melted ice thus increases, resulting in more pollutants adhering to the soil particles which dissolved into the river water [7][8][9]. For instance, in soil aggregates, the opening of the clay lattice will release more fixed ammonia-nitrogen (NH 4 + -N) [10]. (2) The solubilities of the various ions and the factors that affect solubility are distinctive, which leads to large variations The export processes of water and pollutants from the soil into streams during the thawing period are illustrated in Figure 1. The catchment was divided into hydrological units. In each hydrological unit, the melted ice in the soil was modeled to move laterally from soil to stream along a continuous channel formed by the broken soil aggregates above the frozen layer.
The lateral flow flux, q, in the soil above the frozen layer at position x was calculated by the following equation: where h is the height of the migrating water layer (the distance between the liquid surface and the frozen layer), K is the effective hydraulic conductivity of the soil (i.e., the conductivity of the continuous channel formed by the broken soil aggregates during the thawing process), i is the topographic slope, and w is the width of the flow surface at position x. Equation (2) was used to describe the change in the width along the slope: where w c and w b are the parameters used to describe the gradual convergence of the flow along the slope. The lateral flow flux, q, in the soil above the frozen layer at position x was calculated by the following equation: where ℎ _ is the height of the migrating water layer (the distance between the liquid surface and the frozen layer), K is the effective hydraulic conductivity of the soil (i.e., the conductivity of the continuous channel formed by the broken soil aggregates during the thawing process), i is the topographic slope, and w is the width of the flow surface at position x. Equation (2) was used to describe the change in the width along the slope: where wc and wb are the parameters used to describe the gradual convergence of the flow where f is the ratio of the mobile water to the total water. According to the law of the conservation of mass, the change of water storage at position x was calculated as follows: where R is the source sink term, which is the water discharged into the thawed soil from the upper boundary (i.e., rainfall) in the hydrological unit, and S is the water storage capacity of the hydrological unit. Equation (3) was substituted into Equation (4) to obtain Equation (5): The quadratic term in Equation (5) was then expanded to obtain: As the changes in w with respect to x are small, ∂w ∂x was ignored, meaning that Equation (6) was transformed into , Equation (7) was expressed as: Here, S was linearly approximated as: where w and S are the average width and water storage capacity of the liquid flow area, respectively; p is between zero and one; and D is the thickness of the soil layer above the frozen layer. Equation (9) was incorporated into Equation (8) to obtain: The soil water movement equation described by the water storage S, was obtained through simplification: The upper boundary is the recharge boundary. During the soil thawing period, most of the rainfall infiltrated into the soil, as the rainfall amounted to <10 mm in most of the rainfall events and runoff did not occur. The lower boundary is the impervious layer formed by the un-melted ice layer. The water storage capacities at x = L and x = 0 are calculated as follows: S(0, t) = 0 for x = 0 (13)

Simulation of the Pollutant Migration and Transformation in Thawing Soil
During the soil thawing period, melting ice was the main source of movable water in the thawing layers, while the pollutants that adhered to the soil particles dissolved in the water. The concentration of the pollutants in the moving water was affected by the mass of the dissolved pollutants, the liquid water content and the soil flow rate. In order to grasp the mechanism of the pollutant migration and the interactions between the soil and liquid water, the migration and transformation processes from the soil into streams during the thawing period were described using three methods independently.
(1) Method 1 In the thawed soil, the pollutant concentration in the moving water was mainly affected by the pollutant mass which dissolved in to the moving water. The pollutant concentration in the thawed soil was expressed as follows: where c c and c m are the pollutant concentrations in the mobile water and soil, respectively, and D s (S) is the hydrodynamic dispersion coefficient of the pollutants. Because pollutants mainly move in free liquid water, D s (S) was given a fixed value. q is the soil liquid water flux, which was calculated using Equation (3), and K m and α are the parameters characterizing the movement of pollutants from soil to mobile water.
(2) Method 2 The processes whereby the pollution adhering to soil particles dissolved in water and the transformation process in liquid water were described using lumped first-order kinetic equations. These two processes formed the source/sink term in the solute transport equation: where k ac = 10 −ΩS k ms is the kinetic parameter modified by S and Ω, and Ω characterizes the bending degree of the soil. k de is the first-order kinetic coefficient; it characterizes the changes in the pollutant concentrations in moving water.
(3) Method 3 The pollutant dissolving process from soil particles to water was described as follows: The pollutant transport and transformation processes were simulated by: where β k , η k , and k sk are the kinetic coefficients.

Simulation of the Hydrology, Migration and Transformation Processes of Pollutants in Streams
The water flow at the outlet of the catchments originates from exported water in soil. The Manning formula was used to simulate the flow rate: where A is the area of the flow section, C is the Chezy coefficient, R is the hydraulic radius, i is the slope of the stream, and n is the roughness of the channel. The transformation process of the pollutant in the streams was described using the first-order dynamic formula: where C (out) is the concentration of pollutants at the outlet of the catchment, c 0 is the pollutant concentration in the water exported from the soil, k st is the first-order kinetic coefficient of the pollutant transformation in the stream, l is the length of the pollutant migration path (half of the length of the stream in the catchment), and u is the flow velocity, which was calculated using Equation (18).

Field Experiment
In order to calibrate and verify the parameters of the model, field experiments were conducted in the Heidingzi river basin of Shuangyang District, Changchun City, Jilin Province, during the soil thawing period from 2016 to 2018. This basin covers 74.4 km 2 , with an altitude range of 114-233 m. Paddy and corn fields account for 14.5% and 67.2% of the total river basin area, respectively, while forests account for 11.2%. The basic information of the experimental catchments is shown in Table 1 and Figure 2. The monitoring data show that the groundwater level was lower than the bottom of the river during the freeze-thaw periods. During the thawing period, the water and pollutants in the river were all exported from the soil, and the rural domestic pollutants in the catchment were mainly discharged into farmland, rather than being directly discharged into the river. The maximum frozen soil depth in the experimental area was 158 cm and the minimum temperature during the freezing period (January) was −38.4 • C. The field experiments were conducted during the thawing periods of the frozen soil (i.e., from after the snow had completely melted until the paddy fields started to flood) in 2016, 2017 and 2018. Monitoring sections were set up at the outlets of three catchments dominated by corn fields, one catchment dominated by paddy fields, and at the outlet of the basin (Figure 2). A Doppler flow meter was used to measure the flow rates. Water samples were taken at these five locations every day during the soil thawing period.
During the thawing period, soil sampling was conducted every 5-7 days. Six locations were selected in each of the above five areas. At each location, soil samples were taken from 0 to 160 cm depth at 10 cm intervals. The frost front depths were estimated based on the soil rigidity. After measuring the water content, the soil samples were air dried, mixed with distilled water at a ratio of 1:5, and shaken for 6 h. The filtered solutions and water samples taken from the streams were used to measure the NH 4 + -N, nitrate-nitrogen (NO 3 -N), and soluble phosphorus (SP) concentrations using an automatic chemical analyzer (Cleverchem 200, DeChem-Tech Company, Hamburg City, German).
tions were selected in each of the above five areas. At each location, soil samples were taken from 0 to 160 cm depth at 10 cm intervals. The frost front depths were estimated based on the soil rigidity. After measuring the water content, the soil samples were air dried, mixed with distilled water at a ratio of 1:5, and shaken for 6 h. The filtered solutions and water samples taken from the streams were used to measure the NH4 + -N, nitratenitrogen (NO3 -N), and soluble phosphorus (SP) concentrations using an automatic chemical analyzer (Cleverchem 200, DeChem-Tech Company, Hamburg City, German).

Parameter Calibration and Model Evaluation
The model mainly simulated the processes by which water and pollutants were exported from soil into streams. The measured values of the soil pollution concentrations, frost front depths and rainfall during the thawing process were used as the input data. The atmospheric temperature and rainfall during the thawing period of the frozen soil in 2016, 2017 and 2018 are shown in Figure 3.

Parameter Calibration and Model Evaluation
The model mainly simulated the processes by which water and pollutants were exported from soil into streams. The measured values of the soil pollution concentrations, frost front depths and rainfall during the thawing process were used as the input data. The atmospheric temperature and rainfall during the thawing period of the frozen soil in 2016, 2017 and 2018 are shown in Figure 3.
tions were selected in each of the above five areas. At each location, soil samples were taken from 0 to 160 cm depth at 10 cm intervals. The frost front depths were estimated based on the soil rigidity. After measuring the water content, the soil samples were air dried, mixed with distilled water at a ratio of 1:5, and shaken for 6 h. The filtered solutions and water samples taken from the streams were used to measure the NH4 + -N, nitratenitrogen (NO3 -N), and soluble phosphorus (SP) concentrations using an automatic chemical analyzer (Cleverchem 200, DeChem-Tech Company, Hamburg City, German).

Parameter Calibration and Model Evaluation
The model mainly simulated the processes by which water and pollutants were exported from soil into streams. The measured values of the soil pollution concentrations, frost front depths and rainfall during the thawing process were used as the input data. The atmospheric temperature and rainfall during the thawing period of the frozen soil in 2016, 2017 and 2018 are shown in Figure 3. The model parameters were determined using the following objective function: where m q is the number of monitored variables (i.e., the flow, NH 4 + -N, NO 3 -N and SP concentrations at the outlets of the catchments), n qj is the monitoring number of the j th variable, and O j (x, t i ) is the monitoring value of the j th variable at monitoring position x, time t i . P j (x, t i , b) is the simulation result with the parameter matrix b, and v j is the weight of variable J, which was used to reduce the influence of the magnitude and number of different monitored variables on the objective function. The Nash Sutcliffe coefficient, N; the relative total error, F B ; and the information entropy, E, were used to evaluate the model performance. N was used to evaluate the goodness-of-fit of the model: where O t and P t are the observed and simulated values at time t, respectively, O is the average of the observed values, and n is the number of observations. F B was used to evaluate the systematic deviation between the simulated and observed values: Here, α = P O , where P and O are the averages of the observed and simulated values, respectively. Equation (24) was expressed as follows: Here, F B = 0 when α = 1, which indicates that there was no systematic error. As information entropy can characterize a time series, this was used to evaluate the interactive impacts of the soil hydrological and transport processes on the pollutant concentrations. The entropy was calculated as follows: where p i is the probability of pollutant concentrations at concentration level i. For each pollutant during each soil thawing period, the concentrations in the stream were divided into 20 levels (i.e., n = 20). A larger E indicated a more disposed distribution of concentrations.

Evaluation of the Simulation of Soil Flow Processes
The model parameters were calibrated using the observations taken during the soil thawing periods in 2016 and 2017. The calibrated parameters are shown in Table 2. The flow and transport processes in the soil and streams during the thawing period in 2018 were simulated using the calibrated parameters to evaluate the model performances. Figure 4 shows a comparison of the simulated and observed flow rates and pollutant concentrations at the outlets of catchments dominated by paddy fields and corn fields during the soil thawing periods in 2016, 2017 and 2018. As shown in Figure 3, the average temperatures during the soil thawing periods in 2016, 2017 and 2018 were 7.5, 7.2 and 7.1 • C, respectively. The temperature during the thawing period in 2016, especially at the beginning, was significantly higher than those in 2017 and 2018. Both the simulated and observed water exports were significantly higher in 2016 than those in 2017 and 2018. Compared with 2016, the amounts of water exported from the soils into streams in 2017 and 2018 were 42.9% and 35.2% lower in the paddy-dominated catchment, respectively. Furthermore, the impact of meteorology on the water export from the corn-dominated catchments was significantly higher than that of the paddy-dominated catchment. In the corn-dominated catchments, the water exports in 2017 and 2018 were 50.7% and 51.8% lower than in 2016, respectively.   As shown in Tables 1 and 2, the calibrated effective hydraulic conductivity, K, which describes the permeability of the soil, was 1.64 × 10 2 times larger than the hydraulic conductivity in unfrozen soil. The freezing and thawing processes affected the water conductivity of the soil by changing the aggregate structure. If the soil flow process was dominated by the soil's hydraulic conductivity, in the paddy-dominated catchment the flow peaked under the saturated conditions. This was inconsistent with the monitored water  Tables 1 and 2, the calibrated effective hydraulic conductivity, K, which describes the permeability of the soil, was 1.64 × 10 2 times larger than the hydraulic conductivity in unfrozen soil. The freezing and thawing processes affected the water conductivity of the soil by changing the aggregate structure. If the soil flow process was dominated by the soil's hydraulic conductivity, in the paddy-dominated catchment the flow peaked under the saturated conditions. This was inconsistent with the monitored water export process, which indicated that the export of water primarily depended on the amount of melted water. For the corn-dominated catchments, the soil was unsaturated, and the amount of exported water was affected by the initial water storage. The inter-annual variation of the exported water in these areas was more obvious.

As shown in
The simulation effect of the model was better when there was more exported water. This can be attributed to the fact that the model used changes in the soil water storage to describe the process of soil water moving in the thawed layer from the frozen soil. Therefore, compared with the measured values, the model performed best in both cornand paddy-dominated catchments in 2016. In the same period, the paddy-dominated catchment benefited from a larger volume of exported water, and the simulation effect was better than that of corn, with a smaller relative error. Overall, Equation (11) accurately described the process of the water export during the soil thawing period.

Simulation of Pollutant Transport and Transformations in the Paddy-and Corn-Dominated Catchments
The exported flow per unit area in the corn-dominated catchments was 68.4% of that in the paddy-dominated catchment. The lower water flux (velocity) allowed water to stay longer in the soil, meaning that more pollutants adhering to the soil particles dissolved in the water. Therefore, the NH 4 + -N, NO 3 -N and SP concentrations at the outlet of the corn-dominated catchment were higher than those of the paddy-dominated catchment by 12.4%, 22.8% and 14.7%, respectively.
For both the paddy-and corn-dominated catchments, the information entropy was smallest for the SP concentrations. The difference in SP entropy between the paddy-and corn-dominated catchments was highest among the three pollutants. The distribution of SP in the streams was relatively concentrated compared to the other pollutants, whereas NH 4 + -N and NO 3 − -N were more dispersed. This was consistent with the parameter calibration results of the three pollutants ( Table 2). The kinetic coefficients (i.e., k m , k ms , and k sk in Equations 14-16) used to describe the changes in the SP concentrations in the soil were the smallest among the three pollutants. In the paddy-and corn-dominated catchments, the differences in the information entropy of NO 3 − -N and SP were 0.267 and 0.395, respectively. The differences between NH 4 + -N and SP were 0.524 and 0.467, respectively. The difference between the information entropy of NO 3 − -N and SP was smaller than that between NH 4 + -N and SP. The migration abilities of NO 3 − -N and SP have been shown to be notably different [1]. The results presented here also indicate that NH 4 + -N was more affected by the transformation than NO 3 − -N in the mobile water. This interpretation was also consistent with the calibration results of the kinetic parameters shown in Table 2. Figure 4 shows a comparison of the simulated and observed pollutant concentrations at the outlets of the paddy-and corn-dominated catchments during the soil thawing periods of 2016, 2017, and 2018. For the three simulation methods of pollutants, method 1 only considered the influence of the pollutant concentration in the soil on the pollutant concentration in the mobile water. Compared with methods 2 and 3, the N-S coefficients of the three pollutants in method 1 were 20.3% and 13.8% smaller on average, respectively, while the relative deviations increased by 21.1% and 29.2% (Table 3). The accuracy of methods 2 and 3 was significantly higher than that of method 1. The information entropy of method 1 during the calibration period (i.e., the thawing period in 2016 and 2017) and the model evaluation period (i.e., the thawing period in 2018) was relatively close, while those of methods 2 and 3 were significantly different. This indicates that the contact between melting ice and soil particles, and the transformation process of pollutants in mobile water, had significant impacts on the pollutants being dissolved into the stream. This finding is consistent with previous research [22].
The pollutants adhered to the soil particles were in full contact with water as the soil pores were fully saturated, and the concentration of the pollutants in the soil had a greater impact on the export process. Method 1 performed better in this period. The contact between migrated water and the pollutants adhered to the soil decreased with the decreasing soil water content. Meanwhile, as the flow rate decreased, the contact time increased. The pollutant concentrations in the mobile water played a key role in the export of pollutants. Methods 2 and 3 considered these factors at this stage and performed well. As discussed above, methods 2 and 3 were more suitable than method 1 in simulating the transport and transformation processes in thawing soil. Method 2 had the highest N-S coefficient, while method 3 had the smallest systematic deviation.

Simulation of Exported Water and Pollutant Transport and Transformation in the Basin
The exported water per unit area of the basin was close to that of the corn-dominated catchments because of the relatively high proportion of corn in the basin. Table 4 shows that the simulated average value of exported water per unit area in the small basin was closer to the measured value, which was better than the catchments. However, as the river process also involved other hydrological processes such as deep seepage, the fitting effect was not as good as that of the catchments (Figure 4). Therefore, the simulation method needs to be further improved. The concentrations of pollutants at the outlet of the basin were affected by the concentrations of pollutants and the amount of exported water in the two catchments. As is consistent with the paddy-and corn-dominated catchments, method 2 also accurately simulated pollutants at the outlet of the basin. Although the exported water per unit area of the paddy-dominated catchment was greater than that of the corn-dominated catchments, the area of the corn-dominated catchments was 4.62 times larger than the paddy-dominated catchment. Therefore, more water was exported from corn-dominated catchments, such that the concentrations of pollutants at the outlet of the basin should have been closer to the actual values for the corn-dominated catchments. Table 4 shows that the concentrations of pollutants in 2017 and 2018 were consistent with this law. However, in 2016, the pollutant concentrations at the outlet of the basin were closer to those in the paddy-dominated catchment where the concentrations were smaller. This may be attributed to the effect of the hydrodynamic characteristics on the reaction's progress, which caused the pollutants to degrade more [29] when the flow rate was high. In addition, compared with the degradation coefficient in summer (0.13-0.26 d −1 ), as calculated in other studies [30,31], the low temperature in spring likely affected the activities of microorganisms [32], which could explain the smaller k st obtained in this study.

Discussion
The two main factors affecting non-point source pollution are the flushing and leaching of water (from irrigation or precipitation), and the field fertilization and accumulated nutrients (source and sink items) [33]. During the soil freeze-thaw period, the source was the accumulated pollutants, and the water came from the melting of snow and ice. When simulating pollutants, the amount of water affected the contact time and degree between water and pollutants attached to the soil, thereby affecting the export process of pollutants, which was the dominant factor affecting non-point source pollution. This result is consistent with previous studies [34]. However, ignoring the impermeability of frozen soil would make the simulation of the water export in the thawing phase lower, which would affect the accuracy of the simulation of the pollutants. At the initial thawing stage, the amount of water exported was large and contained no pollutants, and the concentration of the pollutants in the soil had a greater impact on the export process. Method 1 performed better at this stage. The mass of pollutants dissolved from the soil decreased with the decreasing water export. The influences of the pollutant transport processes on the pollution concentration increased. This was attributed to the fact that when the amount of the exported water was small, the frozen soil layer caused the water that should have been supplied to the unsaturated soil or groundwater to be discharged into the river. At this stage, the pollutants in the soil and mobile water had sufficient time to interact, and the concentration of the soil and mobile water was the dominant factor. Methods 2 and 3 performed better than Method 1. Therefore, the impact of the contact between soil pollutants and mobile water, and the migration of pollutants in the water on the simulation of non-point source pollutions during thawing cannot be ignored. In addition, using a higher degree of freedom to describe the migration of pollutants in soil and flowing water can improve the goodness-of-fit between the simulated and measured values to a certain extent, but it increased the system deviation (Method 2).

Conclusions
During the soil freezing period, roots and vegetation are reduced, and the biological processes of microorganisms are disturbed, causing the accumulation of pollutants in the soil; when the soil is thawing, the pollution is concentratedly exported into the river. This study considered different underlying surfaces and a small river basin as the research objects and constructed a model to simulate the water and pollution export from soil to stream during the thawing period. The field experiment data collected in Shuangyang District, Changchun City, Jilin Province during the soil thawing periods of 2016, 2017 and 2018 were used to calibrate the model's parameters and evaluate its performance. The main conclusions are as follows: • Numerical simulations of hydrological and pollution transport and transformation behaviors at the small river basin-scale were conducted. The proposed water flow motion equation accurately reflected the flow mechanism of melted ice flowing through channels formed by broken soil aggregates. The constructed model effectively described the pollution export from soil to stream during the thawing period.

•
The concentration of the pollutants peaked when the initial amount of exported water was large. The pollutant transport processes influenced the pollution export more significantly after the soil water was significantly reduced.

•
The concentrations of pollutants in the soil, the contact between the soil pollutants and mobile water, and the migration of pollutants in the water all had a significant impact on the export of pollutants. The main controlling factors affecting the pollutions export from soil to stream changed constantly.
It should be pointed out that in order to accurately evaluate the performance of the model, the measured values of the soil pollutant concentrations during the thawing process were used as input data. However, the cumulative amount of each pollutant before the soil begins to thaw was different, and the effects of the soil freezing on the mineralization and denitrification of nitrogen and phosphorus were different. Future studies will simulate the accumulation of pollutants during the freezing process, such that the model can be better coupled with the existing hydrological and pollutant models.