Migration Law of LNAPLs in the Groundwater Level Fluctuation Zone Affected by Freezing and Thawing

: Freezing and thawing can cause dynamic ﬂuctuations of the groundwater level, resulting in the migration and retention of LNAPLs. However, this process is difﬁcult to observe visually, and a suitable simulation method for its quantitative calculation is lacking. In this study, a numerical simulation is established by coupling the HYDRUS-1D software and the TOUGH program to realize dynamic simulation of the entire process of soil temperature changes, water migration, water level ﬂuctuation, and redistribution of LNAPLs during the freeze–thaw process. The results of the study show that the process of soil freezing and thawing causes water migration, which in turn causes groundwater level ﬂuctuation, leading to the migration and redistribution of LNAPLs within the water level ﬂuctuation zone. In this process, the soil particle size and porosity control the response degree and speed of the water level under freezing and thawing and the spatiotemporal distribution of LNAPLs by affecting the soil temperature, capillary force, and water migration. The migration ability of free LNAPLs is determined by their own density and viscosity; the retention of residual LNAPLs is affected by soil porosity and permeability as well as LNAPL viscosity. The results of this study can not only be used to develop a simulation method for the migration and retention mechanism of LNAPLs in cold regions but also serve as a scientiﬁc and theoretical basis for LNAPL pollution control in seasonal frozen soil regions. close to 0.05 are distributed in the descending range. This indicates that a considerable part of the free-state LNAPLs changes into the residual state and remains in the process of the water level drop. In addition, there are some LNAPLs below the water level during the water level drop. This proves that LNAPLs will preferentially displace water when the water level drops [40].


Introduction
The exploration, development, transportation, and utilization of petroleum are often accompanied by serious groundwater and soil pollution problems. The main pollutants are petroleum hydrocarbons in the form of light nonaqueous phase liquids (LNAPLs), which are carcinogenic and teratogenic to the human body and pose a continuous threat to the ecological environment [1]. After the pollution of LNAPLs leaks on the surface, under the action of gravity and capillary force it will gradually diffuse and migrate to the deep layer of the vadose zone and its surroundings, eventually making it difficult to use groundwater resources and soil [2,3].
Seasonal frozen soil is a type of soil in which the surface layer is frozen in winter and thawed in summer. Approximately 30% of the total land area of the northern hemisphere consists of this type of soil, and there are a large amount of oil and oil fields located in these seasonal frozen soil areas which are contributing to the pollution of LNAPLs in the freeze-thaw environment [4,5]. Compared with non-frozen soil areas, seasonal frozen soil areas experience more than 100 days of surface freezing and thawing processes each of research can only focus on the migration of LNAPLs in the water level fluctuation zone under normal temperature environments and is not suitable for seasonal frozen soil regions [42][43][44][45].
Based on the foundation and deficiencies of related research, in this study, HYDRUS-1D, which can simulate the water and heat transfer in the freeze-thaw process [46], is coupled with the TOUGH program, which can simulate the distribution of LNAPLs with water level fluctuation. Through this coupled simulation method, under the conditions of different media and different types of LNAPLs, the migration and retention law of LNAPLs with the water level fluctuation during freeze-thaw periods is comprehensively analyzed, and the factors affecting it are evaluated. This research is expected to provide the theoretical methods and promotional basis for the numerical simulation of LNAPL migration in the water level fluctuation zone of seasonal frozen soil regions as well as serve as a reference for pollution evaluation, management, and the restoration of LNAPLs.

Soil Column Experiment of Water and Heat Transport
The equipment used in the experiment consists of a freeze-thaw cycle box, a refrigeration system, moisture content and temperature probes, a pressure measuring tube, and a water supply system with a constant water head. The experimental setup is shown in Figure 1.
distribution, partitioning, and fate of NAPL in the vadose zone. However, becau TOUGH model has limitations in the simulation of multiphase flow in sub-zero tem ture environments, this type of research can only focus on the migration of LNAPL water level fluctuation zone under normal temperature environments and is not s for seasonal frozen soil regions [42][43][44][45].
Based on the foundation and deficiencies of related research, in this study, HY 1D, which can simulate the water and heat transfer in the freeze-thaw process coupled with the TOUGH program, which can simulate the distribution of LNAPL water level fluctuation. Through this coupled simulation method, under the condit different media and different types of LNAPLs, the migration and retention LNAPLs with the water level fluctuation during freeze-thaw periods is comprehen analyzed, and the factors affecting it are evaluated. This research is expected to p the theoretical methods and promotional basis for the numerical simulation of L migration in the water level fluctuation zone of seasonal frozen soil regions as w serve as a reference for pollution evaluation, management, and the restoration of LN

Soil Column Experiment of Water and Heat Transport
The equipment used in the experiment consists of a freeze-thaw cycle box, a eration system, moisture content and temperature probes, a pressure measuring tub a water supply system with a constant water head. The experimental setup is sho Figure 1. The test soil column is a plexiglass container with a height of 30 cm, and m sand (average particle size 0.5-0.2 mm) was selected as the representative test med The test soil column is a plexiglass container with a height of 30 cm, and medium sand (average particle size 0.5-0.2 mm) was selected as the representative test medium. A water supply pipe and a Markov bottle were used to supply water at a fixed head to the soil column, and the initial water level was set to 16 cm from the top of the soil.
A total of five monitoring holes were set on the wall of the column at an interval of 4 cm to monitor the soil temperature and moisture content in the soil column in real time. Each monitoring hole is equipped with a CR300 soil temperature and humidity sensor. The piezometer reading was periodically taken to obtain the dynamic data of the water level. After the start of the experiment, the piezometer is read every 10 min, and the temperature and water content are automatically monitored and outputted every 2 min through the sensor.
The experimental soil column was placed in the freeze-thaw cycle box, which adjusts the temperature around the soil column to be constant at 2 • C. The temperature control cover on the top of the soil is used to control the freezing and thawing of the soil from the top to the bottom in stages. The minimum temperature in the freezing period is −35 • C, and the temperature in the thawing period is 5 • C. The temperature control process is shown in Figure 2. soil column, and the initial water level was set to 16 cm from the top of the soil A total of five monitoring holes were set on the wall of the column at an in cm to monitor the soil temperature and moisture content in the soil column in Each monitoring hole is equipped with a CR300 soil temperature and humid The piezometer reading was periodically taken to obtain the dynamic data of level. After the start of the experiment, the piezometer is read every 10 min, an perature and water content are automatically monitored and outputted ev through the sensor.
The experimental soil column was placed in the freeze-thaw cycle box, wh the temperature around the soil column to be constant at 2 °C. The temperatu cover on the top of the soil is used to control the freezing and thawing of the so top to the bottom in stages. The minimum temperature in the freezing period and the temperature in the thawing period is 5 °C. The temperature control shown in Figure 2. The resulting data from the soil column experiment will be used to adjust parameters in the subsequent HYDRUS-1D model. Taking the empirical value drothermal parameters in HYDRUS-1D as the initial value [46], the hydraulic p related to the soil hydraulic conductivity in the medium sand obtained by adj fitting are shown in Table 1, and the parameters related to the soil heat condu shown in Table 2. And for the data of the soil column experiment, please refer t in the Supplementary Materials.

Model Method
HYDRUS-1D and TOUGH are coupled to establish a model to simulate th the freeze-thaw cycle on temperature, moisture migration, water level fluctu LNAPL migration and retention. The coupling process mainly considers the c of temperature and water level between the two codes. The simulation range DRUS-1D model is from the top of the soil column (0 cm) to 20 cm below, and groundwater level is located 16 cm below the column surface. The upper boundaries of the water flow are zero flux boundaries, and the upper and lowe ries of heat transfer are temperature boundaries that change with time. Since T The resulting data from the soil column experiment will be used to adjust and fit the parameters in the subsequent HYDRUS-1D model. Taking the empirical value of the hydrothermal parameters in HYDRUS-1D as the initial value [46], the hydraulic parameters related to the soil hydraulic conductivity in the medium sand obtained by adjusting and fitting are shown in Table 1, and the parameters related to the soil heat conductivity are shown in Table 2. And for the data of the soil column experiment, please refer to Table S1 in the Supplementary Materials.

Model Method
HYDRUS-1D and TOUGH are coupled to establish a model to simulate the effect of the freeze-thaw cycle on temperature, moisture migration, water level fluctuation, and LNAPL migration and retention. The coupling process mainly considers the consistency of temperature and water level between the two codes. The simulation range of the HYDRUS-1D model is from the top of the soil column (0 cm) to 20 cm below, and the initial groundwater level is located 16 cm below the column surface. The upper and lower boundaries of the water flow are zero flux boundaries, and the upper and lower boundaries of heat transfer are temperature boundaries that change with time. Since TOUGH is not suitable for simulation at a temperature below 0 • C, the TOUGH simulation range is the area within −10 cm to −30 cm in the soil column experiment, and the soil temperature change within this depth range is consistently above 0 • C. And order to comply with the aforementioned soil column experiment and the HYDRUS-1D numerical model, the top is set as the Dirichlet boundary, with constant atmospheric pressure, no water infiltration, and zero water flux. The bottom is set as a waterproof boundary condition. The initial temperature of the topper and lower are set to 6.55 • C and 7 • C, respectively, according to the experimental results, and the temperature changes with time during the simulation. And since the heat and water migration in the freezing and thawing process are mainly reflected in the vertical direction, the boundary conditions in the coupling model are assumed to be horizontally closed [6][7][8][9][10].The left boundary (x = 0-0.05 m) and the right boundary (x = 0.85-0.90 cm) of the simulation area are both physical boundaries, which are set as zero flux boundary conditions [39].
Coupled simulation is achieved through the transfer of heat and water volume time by time between HYDRUS-1D and TOUGH, and the temperature and water content at each depth are verified by measured data. The interface for the water flux interaction of the coupled model is the layer above the water level fluctuation in the TOUGH model. In order to realize the mutual coupling of water flow, 16 cells in this layer are used as the water injection and pumping cells to control the water level changes, and the changed flux is entered into the source and sink of the cell. Figure 3 shows the water level control obtained by TOUGH. On the other hand, the interaction of heat in the coupled model is achieved jointly by each layer of cells in TOUGH. The heat source-sink items are added to each cell after the division in TOUGH, and the corresponding conversion relationship between the rate of heat addition and the temperature change at each representative layer is calculated. Thus, the coupling of water and heat between the two models is realized. change within this depth range is consistently above 0 °C. And order to comp aforementioned soil column experiment and the HYDRUS-1D numerical mo is set as the Dirichlet boundary, with constant atmospheric pressure, no water and zero water flux. The bottom is set as a waterproof boundary condition. temperature of the topper and lower are set to 6.55 °C and 7 °C, respectively, a the experimental results, and the temperature changes with time during the And since the heat and water migration in the freezing and thawing process reflected in the vertical direction, the boundary conditions in the coupling mo sumed to be horizontally closed [6][7][8][9][10].The left boundary (x = 0-0.05 m) an boundary (x = 0.85-0.90 cm) of the simulation area are both physical bounda are set as zero flux boundary conditions [39].
Coupled simulation is achieved through the transfer of heat and water v by time between HYDRUS-1D and TOUGH, and the temperature and water each depth are verified by measured data. The interface for the water flux in the coupled model is the layer above the water level fluctuation in the TOUGH order to realize the mutual coupling of water flow, 16 cells in this layer are water injection and pumping cells to control the water level changes, and the ch is entered into the source and sink of the cell. Figure 3 shows the water level tained by TOUGH. On the other hand, the interaction of heat in the couple achieved jointly by each layer of cells in TOUGH. The heat source-sink items to each cell after the division in TOUGH, and the corresponding conversion r between the rate of heat addition and the temperature change at each represen is calculated. Thus, the coupling of water and heat between the two models is , θv is the volumetric vapor content expressed as an In the water and heat transfer model, The modified Richards' equation [43] is used as the basic equation of water migration: where θ u is the volumetric unfrozen water content (L 3 L −3 ) (= θ + θ u ), θ is the volumetric liquid water content (L 3 L −3 ), θ v is the volumetric vapor content expressed as an equivalent water content ( , h is the pressure head (L), T is the temperature (K), and S is a sink term (T −1 ) usually accounting for root water uptake. Water flow in Equation (1) is assumed to be caused by five different processes with corresponding hydraulic conductivities. The first three terms on the righthand side of Equation (1) represent liquid flows due to a pressure head gradient (K Lh , [LT −1 ]), gravity, and a temperature gradient (K LT , [L 2 T −1 K −1 ]), respectively. The next two terms represent vapor flows due to pressure head (K vh , [LT −1 ]) and temperature (K vT , [L 2 T −1 K −1 ]) gradients, respectively. The soil hydraulic conduction module uses the Van Genuchten-Mualem equation, which is generally recognized in soil water calculations. The basic equation of soil heat transfer [28] is given as follows.
where the first term on the left-hand side represents changes in the energy content, and the second and third terms represent changes in the latent heat of the frozen and vapor phases, respectively. The terms on the right-hand side represent soil heat flow by conduction, convection of sensible heat with flowing water, transfer of sensible heat by diffusion of water vapor, transfer of latent heat by diffusion of water vapor, and uptake of energy associated with root water uptake, respectively. The phase change between water and ice is controlled by the generalized Clapeyron equation, which defines a relationship between the liquid pressure head and temperature when ice is present in the porous material [32]. Hence, the unfrozen water content can be derived from the liquid pressure head as a function of temperature alone when ice and pure water coexist in the soil [43]. Further, ice volume fraction is determined by the difference between initial water content and liquid water content. The volumetric heat capacity of the soil, (2) is defined as the sum of the volumetric heat capacities of the solid (C n ), liquid (C w ), vapor (C v ), and ice (C i ) phases multiplied by their respective volumetric fractions: Furthermore, in Equation (2), L o is the volumetric latent heat of vaporization of , and L f is the latent heat of freezing (J kg −1 , L 2 T −2 ) (approximately 3.34 × 10 5 J kg −1 ).
In the simulation of LNAPL migration and retention with water level fluctuation, the mass and energy conservation equations [37] are used as the basic governing equations of the components in the multiphase flow (i.e., water, gas, and NAPL phases). Parker's equation [39] is used as the relative permeability coefficient equation in the multiphase flow problem. Van Genuchten's equation [39] is used as the capillary pressure governing equation of the NAPL-water phase. Finally, the migration and flow equation of free LNAPLs in the soil medium [38] is shown in Equation (4).
where φ is the porosity, S N is the NAPL phase saturation, ρ N is the NAPL phase density, and V is the Darcy flow velocity of the immiscible phase. The last term I on the right side of the equation is the source (sink) term (=ρ N Q N + E), Q N is the volume of NAPL existing or leaked in the unit time, and E represents the mass exchange between phases. Such source and sink terms are not considered in this model, so I = 0.

Law of Soil Water and Heat Transfer during Freezing and Thawing
During the freeze-thaw experiment, the water content within 4 cm of the top of the soil column changes drastically, becoming weaker with increasing depth (Figure 4(c1)). This is mainly related to the mutual transformation of the water and ice phases. During the freezing process, as shown in Figure 4(b1,c1), the upper part of the soil body freezes first, decreasing the water content. The negative temperature transfers from top to bottom and the water migrates upward, causing the water level to drop, as shown in Figure 4(d1). In the process of melting, the ice phase on the upper part of the soil melts, the water content increases, and the water moves downward, which increases the groundwater level. The HYDRUS-1D model also shows the changes of temperature, water content and water level in the medium-sand medium during the freezing and thawing process (Figure 4(a'2-d'2)), and its water and heat transfer law is consistent with the results of the soil column experiment: Moisture migration due to freezing and thawing of the upper soil is responsible for changes in the water table.
Water 2022, 14, x FOR PEER REVIEW 7 of 15

Law of Soil Water and Heat Transfer during Freezing and Thawing
During the freeze-thaw experiment, the water content within 4 cm of the top of the soil column changes drastically, becoming weaker with increasing depth (Figure 4(c1)). This is mainly related to the mutual transformation of the water and ice phases. During the freezing process, as shown in Figure 4(b1,c1), the upper part of the soil body freezes first, decreasing the water content. The negative temperature transfers from top to bottom and the water migrates upward, causing the water level to drop, as shown in Figure 4(d1). In the process of melting, the ice phase on the upper part of the soil melts, the water content increases, and the water moves downward, which increases the groundwater level. The HYDRUS-1D model also shows the changes of temperature, water content and water level in the medium-sand medium during the freezing and thawing process (Figure 4(a'2d'2)), and its water and heat transfer law is consistent with the results of the soil column experiment: Moisture migration due to freezing and thawing of the upper soil is responsible for changes in the water table. As shown in Figure 4(1,2), In the HYDRUS-1D model, the changes of soil temperature and water content at each depth are similar to the soil column experiment results. Although the water level changes stepwise, the overall trend and range of change are consistent with the soil column experiment results. In order to evaluate the accuracy of the Hydrus-1D model, two indicators, the percentage of deviation (PBIAS) and the root mean square error (RMSE), were used to quantitatively evaluate the simulation data of HY-DRUS-1D for medium-sand media. The PBIAS absolute value of the water level and temperature data in the HYDRUS-1D model is less than 10%, and RMSE is within 0.6, indicating that the water level and soil temperature simulation results are good, and they are As shown in Figure 4(1,2), In the HYDRUS-1D model, the changes of soil temperature and water content at each depth are similar to the soil column experiment results. Although the water level changes stepwise, the overall trend and range of change are consistent with the soil column experiment results. In order to evaluate the accuracy of the Hydrus-1D model, two indicators, the percentage of deviation (PBIAS) and the root mean square error (RMSE), were used to quantitatively evaluate the simulation data of HYDRUS-1D for medium-sand media. The PBIAS absolute value of the water level and temperature data in the HYDRUS-1D model is less than 10%, and RMSE is within 0.6, indicating that the water level and soil temperature simulation results are good, and they are overall consistent with the experimental results and have high accuracy. Figure 5a-f show the changes in the water level and temperature obtained by the HYDRUS-1D model with different media. The response rate of soil temperature conduction and water level to the freeze-thaw process is in the order of silty clay > fine sand > medium sand > coarse sand. Compared with the other three media, the soil temperature at the middle depth level of silty clay responds more quickly to the temperature changes in the upper part of the soil. This is because the soil solid particles are arranged closer together, allowing the upper temperature to diffuse downward during the freeze-thaw process [16] and thus improving the overall heat transfer properties of the soil. Moreover, because the temperature of silty clay soil changes rapidly, the capillary force of the upper part of the soil also changes faster, and the time of drop and rise of the water level is earlier than that in the case of coarse sand, medium sand, and fine sand. overall consistent with the experimental results and have high accuracy. Figure 5a-f show the changes in the water level and temperature obtained by the HYDRUS-1D model with different media. The response rate of soil temperature conduction and water level to the freeze-thaw process is in the order of silty clay > fine sand > medium sand > coarse sand. Compared with the other three media, the soil temperature at the middle depth level of silty clay responds more quickly to the temperature changes in the upper part of the soil. This is because the soil solid particles are arranged closer together, allowing the upper temperature to diffuse downward during the freeze-thaw process [16] and thus improving the overall heat transfer properties of the soil. Moreover, because the temperature of silty clay soil changes rapidly, the capillary force of the upper part of the soil also changes faster, and the time of drop and rise of the water level is earlier than that in the case of coarse sand, medium sand, and fine sand. The magnitude of the water level change in different media is in the order of fine sand > silty clay > medium sand > coarse sand. For the same temperature gradient, the smaller the particle size of the soil, the better is its water holding capacity and the greater is the capillary force produced during freezing. Therefore, in the freezing process, the amount of water migration on the soil profile in the fine sand medium is greater than that in the medium sand and coarse sand, which increases the water level change.
Compared with fine sand, silty clay has smaller particles; however, owing to its lower residual moisture content and greater soil porosity, silty clay has lower capillary force generated during the freeze-thaw process. Consequently, the water migration and the water level fluctuation range on the soil profile are smaller than those in the case of fine sand. In addition, the water level change range in silty clay is larger than that in coarse sand and medium sand. The reason is similar to that for fine sand: the smaller the soil particle size, the better is the water holding capacity and the greater is the capillary force generated during freezing.
Taken together, in the same temperature gradient, the capillary force generated by the medium during the freeze-thaw process and the thermal conductivity of the medium are important factors that affect the water migration, range, and speed of the water level change. Figure 6 shows the migration and retention of LNAPLs with water level fluctuation (TOUGH model with toluene as the LNAPL and medium sand as the medium). The residual saturation of the LNAPL phase is 0.05; therefore, when the saturation is less than 0.05, the LNAPLs will remain under the control of capillary force and transform into a residual state. As the water level drops, the free-state LNAPLs migrate downward; simultaneously, their saturation continues to decrease and gradually expands to the entire level in the horizontal direction. When the LNAPLs drop to the lowest position with the water level, a part of the LNAPLs with saturation close to 0.05 are distributed in the descending range. This indicates that a considerable part of the free-state LNAPLs changes into the residual state and remains in the process of the water level drop. In addition, there are some LNAPLs below the water level during the water level drop. This proves that LNAPLs will preferentially displace water when the water level drops [40]. temperature change obtained by the silty clay model. The magnitude of the water level change in different media is in the order of fine sand > silty clay > medium sand > coarse sand. For the same temperature gradient, the smaller the particle size of the soil, the better is its water holding capacity and the greater is the capillary force produced during freezing. Therefore, in the freezing process, the amount of water migration on the soil profile in the fine sand medium is greater than that in the medium sand and coarse sand, which increases the water level change.

Migration and Retention Process of LNAPLs Affected by Freezing and Thawing
Compared with fine sand, silty clay has smaller particles; however, owing to its lower residual moisture content and greater soil porosity, silty clay has lower capillary force generated during the freeze-thaw process. Consequently, the water migration and the water level fluctuation range on the soil profile are smaller than those in the case of fine sand. In addition, the water level change range in silty clay is larger than that in coarse sand and medium sand. The reason is similar to that for fine sand: the smaller the soil particle size, the better is the water holding capacity and the greater is the capillary force generated during freezing.
Taken together, in the same temperature gradient, the capillary force generated by the medium during the freeze-thaw process and the thermal conductivity of the medium are important factors that affect the water migration, range, and speed of the water level change. Figure 6 shows the migration and retention of LNAPLs with water level fluctuation (TOUGH model with toluene as the LNAPL and medium sand as the medium). The residual saturation of the LNAPL phase is 0.05; therefore, when the saturation is less than 0.05, the LNAPLs will remain under the control of capillary force and transform into a residual state. As the water level drops, the free-state LNAPLs migrate downward; simultaneously, their saturation continues to decrease and gradually expands to the entire level in the horizontal direction. When the LNAPLs drop to the lowest position with the water level, a part of the LNAPLs with saturation close to 0.05 are distributed in the descending range. This indicates that a considerable part of the free-state LNAPLs changes into the residual state and remains in the process of the water level drop. In addition, there are some LNAPLs below the water level during the water level drop. This proves that LNAPLs will preferentially displace water when the water level drops [40]. When the water level rises, the free LNAPLs also rise with the water level, and, at the same time, their saturation continues to decrease. In this process, residual LNAPLs con tinuously stay below the groundwater surface. Therefore, the drop and rise of the wate level during the entire freeze-thaw process will cause the original LNAPLs to undergo a significant redistribution process. Furthermore, the pollution range of LNAPLs continue to expand both vertically and horizontally with decreasing water level, while the pollu tion range of the area above the initial water level increases as the water level increases During the whole process, the vertical fluctuation range of the water level is approxi mately 3.3 cm, and the vertical pollution range caused by the migration and retention o LNAPLs is approximately 6 cm.

Migration and Retention Process of LNAPLs Affected by Freezing and Thawing
Taken together, the mechanism of the freeze-thaw process on the migration and re tention of LNAPLs is as follows: soil temperature changes cause water phase changes resulting in water migration, which in turn causes water level fluctuations and the migra tion and retention of LNAPLs within the water level fluctuation zone. The temperature a the top of the soil is the most fundamental controlling factor in this process. During th freeze-thaw cycle, when LNAPLs rose and fell with the water level, a considerable par of the LNAPLs was retained and transformed into a residual state. Therefore, the effect o retention cannot be ignored. Both migration and retention effects lead to the continuou expansion of the pollution range of LNAPLs in the water level fluctuation zone, resulting in the redistribution of LNAPLs with the water level. Therefore, compared with that in areas not affected by freeze-thaw, the pollution range of LNAPLs in the water level fluc tuation zone will be larger in the seasonal frozen soil area.

Analysis of Factors Affecting the Migration of LNAPLs
After successfully simulating the migration and retention model of LNAPLs in the freeze-thaw water level fluctuation zone by combining HYDRUS-1D and TOUGH, th model is extended to different media conditions and different types of LNAPLs to study their effect on the migration and retention of LNAPLs. Table 3 shows the comparison o model results with different media conditions and different types of LNAPLs.  When the water level rises, the free LNAPLs also rise with the water level, and, at the same time, their saturation continues to decrease. In this process, residual LNAPLs continuously stay below the groundwater surface. Therefore, the drop and rise of the water level during the entire freeze-thaw process will cause the original LNAPLs to undergo a significant redistribution process. Furthermore, the pollution range of LNAPLs continues to expand both vertically and horizontally with decreasing water level, while the pollution range of the area above the initial water level increases as the water level increases. During the whole process, the vertical fluctuation range of the water level is approximately 3.3 cm, and the vertical pollution range caused by the migration and retention of LNAPLs is approximately 6 cm.
Taken together, the mechanism of the freeze-thaw process on the migration and retention of LNAPLs is as follows: soil temperature changes cause water phase changes, resulting in water migration, which in turn causes water level fluctuations and the migration and retention of LNAPLs within the water level fluctuation zone. The temperature at the top of the soil is the most fundamental controlling factor in this process. During the freeze-thaw cycle, when LNAPLs rose and fell with the water level, a considerable part of the LNAPLs was retained and transformed into a residual state. Therefore, the effect of retention cannot be ignored. Both migration and retention effects lead to the continuous expansion of the pollution range of LNAPLs in the water level fluctuation zone, resulting in the redistribution of LNAPLs with the water level. Therefore, compared with that in areas not affected by freeze-thaw, the pollution range of LNAPLs in the water level fluctuation zone will be larger in the seasonal frozen soil area.

Analysis of Factors Affecting the Migration of LNAPLs
After successfully simulating the migration and retention model of LNAPLs in the freeze-thaw water level fluctuation zone by combining HYDRUS-1D and TOUGH, the model is extended to different media conditions and different types of LNAPLs to study their effect on the migration and retention of LNAPLs. Table 3 shows the comparison of model results with different media conditions and different types of LNAPLs.
The response of water level and LNAPL position to temperature changes during freeze-thaw is affected by the medium condition; the response speed in different media is in the order of fine sand < silty clay < medium sand < coarse sand. This indicates that the smaller the soil particles, the smaller the porosity. That is, with closer arrangement of the soil particles, the heat conduction properties of the soil improves and the hysteresis of the soil temperature change at the middle depth of the soil decreases. As a result, water migration can quickly respond to changes in the top surface temperature of the soil during the freeze-thaw process, and the start of the fluctuation of free LNAPLs near the water level occurs earlier. The water level fluctuation and the final pollution depth range of LNAPLs in different media are in the order of fine sand > silty clay > medium sand > coarse sand. This phenomenon should be related to the porosity of the soil. The smaller the porosity of the soil, the greater is the capillary force generated in the unsaturated zone during freezing, resulting in the increase of the amount of water migration on the soil profile as well as the ranges of water level fluctuation and LNAPL pollution.
The rate of LNAPL migration and expansion is related to the porosity, permeability, and properties of the medium. With higher porosity and permeability of the soil or lower viscosity of LNAPLs, the migration speed of LNAPLs as well as the distribution thickness of free-state LNAPLs that can move with the water level increase. Conversely, with lower soil porosity and permeability or greater viscosity of LNAPLs, the migration speed of LNAPLs in the vertical direction decreases and more LNAPLs will expand and be retained in the lateral direction. Thus, the maximum saturation of LNAPLs becomes lower at the end of the freeze-thaw process. Therefore, at the lowest water level, in the fine sand medium with the least porosity, the residual LNAPLs near the initial position have the largest retention range. Among different LNAPLs, there is less residual chloroethane in the initial position because of the lower viscosity of chloroethane.
In addition, some LNAPLs exist below the water level during the freezing period, and from the perspective of saturation, this part of the LNAPLs is not in a stagnant state. This is possibly because LNAPLs will preferentially displace water when the water level decreases [42]. In each medium, this part of LNAPLs does not have much difference in the range below the water level; however, comparing the three LNAPLs, the denser chloroethane has a wider range below the water level. Therefore, LNAPLs with high density can more easily displace water below the water level.
In the same TOUGH model, each medium is simulated at a constant temperature, only changing the non-isothermal mode to the isothermal mode. The soil temperature was constant at 22 • C, with no change in other conditions such as the water level. The results of isothermal simulation are compared with the results during the freeze-thaw process. The comparison is shown in Table 4. The migration distribution and pollution range of LNAPLs in both isothermal and non-isothermal models do not show any considerable difference. The only difference is that the maximum saturation of LNAPLs in the isothermal condition is approximately 0.01-0.02 larger than that during the freeze-thaw process; this difference is relatively small. This confirms that temperature does not directly control the migration and retention of LANPLs, but indirectly leads to the redistribution of LNAPLs by affecting water level changes.

Conclusions
In this study, the HYDRUS-1D software and TOUGH program were combined to successfully establish a migration and retention model of LNAPLs with water level fluctuation under the effect of the freeze-thaw cycle. This model can thus be applied to study the interaction mechanism of freezing and thawing temperature-dependent LNAPL redistribution in the water level fluctuation zone. Based on this model, the migration models under different media, LNAPLs, and temperature conditions were established. The following conclusions can be drawn: (1) The temperature changes during freezing and thawing is the most fundamental factor controlling the migration of LNAPLs in the water level fluctuation zone. The LNAPLs near the water level during the freeze-thaw process migrate or remain with the fluctuation of the water level. Thus, the pollution range and diffusion degree of LNAPLs in the seasonal frozen soil area are greater than those in the non-frozen soil area. (2) The particle size and porosity of the medium as well as LNAPL composition type are the main factors that affect the migration of LNAPLs. In general, the finer and more closely packed the medium particles, the faster is the response of soil temperature and water level changes to freezing and thawing. The smaller the porosity, the greater is the capillary force generated above the water level during the freeze-thaw process and the greater are the amount of water migration, the amplitude of the water level, and the migration and diffusion range of LNAPLs. By contrast, the composition type of LNAPLs determines their retention process and phase change. (3) The coupling model in this study can provide a simulation method for the analysis of the migration and retention of LNAPLs in the water level fluctuation zone of seasonal frozen soil regions. The research results on the migration and retention mechanism and the factors affecting LNAPLs can be used as a theoretical basis for LNAPL pollution evaluation and treatment in areas subjected to freeze-thaw cycles.
Considering that T2VOC lacks simulation functions for chemical adsorption and microbial degradation processes, the effect of temperature on chemical and biological processes is not reflected and requires further research.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/w14081289/s1, Table S1: Result data of water and heat transfer in medium sand medium soil column experiment, Table S2: HYDRUS-1D simulation results.