A New Mechanism of Canopy Effect in Unsaturated Freezing Soils

Canopy effect refers to the phenomenon where moisture accumulates underneath an impervious cover. Field observation reveals that canopy effect can take place in relatively dry soils where the groundwater table is deep and can lead to full saturation of the soil immediately underneath the impervious cover. On the other hand, numerical analysis based on existing theories of heat and mass transfer in unsaturated soils can only reproduce a minor amount of moisture accumulation due to an impervious cover, particularly when the groundwater table is relatively deep. In attempt to explain the observed canopy effect in field, this paper proposes a new mechanism of moisture accumulation in unsaturated freezing soils: vapour transfer in such a soil is accelerated by the process of vapour-ice desublimation. A new approach for modelling moisture and heat movements is proposed, in which the phase change of evaporation, condensation and de-sublimation of vapor flow are taken into account. The computed results show that the proposed model can indeed reproduce the unusual moisture accumulation observed in relatively dry soils. The results also demonstrate that soil freezing fed by vapour transfer can result in a water content close to full saturation. Since vapour transfer is seldom considered in geotechnical design, the canopy effect deserves more attention during construction and earth works in cold and arid regions.


Introduction
Significant frost heave damage was reported to airfield runways in an airport in north-western China, with runways broke-down or even shattered at some positions.It was found that a weak zone of subgrade soil formed along the centreline of the runway where the water content was observed to be close to full saturation (Li et  al., 2014 [1] ).The airport is located in cold and arid northwestern China.The typical annual precipitation in this region is as low as 328 mm, and the annual evaporation can be as high as 1438 mm.Moisture accumulation due to rainwater infiltration was considered unlikely as the pavement was much intact and impervious during the previous raining season (summer).The annual minimum and maximum temperature are around -20 ºC and 35 ºC, respectively.The groundwater table remains more than 20 m deep, eliminating the possibility of capillary action influence on the water content of the surface soil.The question is then: where is the water in the subgrade soil coming from?
According to Li et al (2014) [1] , the pavement acts as a 'pot cover' and prevents evaporation, leading to moisture accumulation under the pavement.They used the terminology of 'pot cover effect', to reflect the fact that the moisture migration was driven by a temperature gradient in the ground and prevented by an impervious cover.However, as will be shown in the paper, this terminology does not represent what truly occurred under the airfield runway.The temperature gradient alone cannot be the only driving force for the moisture accumulation.There must exist some other mechanism that enables the moisture to accumulate, most likely in form of vapour transfer when the groundwater table is 20m deep.The biological terminology of canopy effect is borrowed in this paper, to reflect the dual mechanisms of vapour transfer.Canopy effect in plant biology refers to the phenomenon where crowns of individual trees or forest act as a special layer of absorption for energy from sunlight and water from tree roots.Significant water is consequently stored in the tree leaves for photosynthesis.The paper will provide a more detailed justification for the use of the terminology to refer to the large amount of moisture accumulation under impervious covers.
The objective of this work is to demonstrate and reproduce the unusual canopy effect observed in relatively cold and dry soils.Numerical simulation of coupled heat and mass transfer based on existing theories is first carried out, in attempt to study the effects of impervious covers on the moisture distribution in unsaturated soils.The limitations of existing theories are demonstrated through numerical examples.A new mechanism of vapour transfer in unsaturated freezing soils is then proposed, in which the phase changes of evaporation, condensation and desublimation during coupled heat and mass transfer are taken into account.The new model is then used to analyse the canopy effect observed in the airfield runways.Finally, some conclusions are drawn based on the numerical results and discussions.
2 Analysis of canopy effect using the existing theory

Canopy effect under impervious covers
When the ground surface is covered by an impervious cover, as shown in Figure 1, the moisture content profiles under the cover (profile A) and outside the cover (profile B) will be different.In addition, two possible moisture profiles may exist under the impervious cover: profile A 1 for a relatively high groundwater table where capillary action is the dominant mechanism of moisture transfer, and profile A 2 for a low groundwater table where vapour transfer is the dominant mechanism.
The occurrence of profile A 1 is straightforward, and indeed the only difference between profile A 1 and profile B is that the moisture content of the former is somewhat higher than the latter, due to the prevention of evaporation of the cover.The occurrence of profile A 2 is however less straightforward and is seldom mentioned in the literature.The case reported by Li et al (2014) [1] refers to profile A 2 in Figure 1.Eigenbrod and Kennepohl (1996) [2] also reported that significant moisture accumulation due to vapour transfer occurred underneath a road pavement and led to significant frost heave damage to the pavement.As profile A 2 represents a deep groundwater table and a dry foundation soil, issues associated with high moisture content are more likely to be ignored in geotechnical design.Therefore, the situation of profile A 2 represents a real threat to engineering practice, and is referred to as the canopy effect in the paper.The main objective here is to reveal the true physical mechanism of profile A 2 and to reproduce it via numerical modelling and laboratory experiments.This paper focuses on the numerical modelling, while the experimental verification will be discussed in a separate paper.

Mathematical model
The theory of coupled moisture and heat transfer in unsaturated soils was pioneered by Philip and de Vries (1957) [3] , who developed a mathematical model to describe the coupled movement of liquid water and water vapor under non-isothermal conditions (PDV model).Milly (1984)  [4] adopted the matric head instead of the total volumetric water content in the PDV model.Nassar and Horton (1989) [5] further extended Milly's work by including osmotic effects on liquid and water vapor movement.Sakai et al. (2009)  [6] presented a coupled model for water and vapor movement which considers condensation and evaporation effects under nonisothermal and low water content conditions.
According to the PDV model, the flow of liquid water and vapor in an unsaturated soil are driven either by a temperature gradient or by a matric head gradient.The governing equation can thus be expressed as: where θ is total water content (m 3 m -3 ); θ L and θ v are liquid water content and equivalent vapor content, respectively (m 3 m -3 ); h is water heat (m); T is temperature (K); z is the spatial coordinate positive upward (m); K Lh (m s -1 ) and K LT (m 2 K -1 s -1 ) are the isothermal and thermal hydraulic conductivities for liquid water fluxes due to gradients in h and T, respectively; K vh (m s -1 ) and K vT (m 2 K -1 s -1 ) are the isothermal and thermal vapor hydraulic conductivities, respectively (Teng et al., 2015 [7] ).
Considering vaporization and condensation processes in unsaturated soil, the governing equation of heat conservation is given as (de Veries, 1958 [8] ; Nassar and Horton, 1992 [5] ): where C p is defined as the sum of the volumetric heat capacity of the unsaturated soil (MJ m −3 K −1 ) and is equal to the volumetric average of the solid (C n = 1.92 MJ m −3 K −1 ), liquid (C L = 4.18 MJ m −3 K −1 ), and air (C v =6.3 KJ m −3 K −1 ) phases (de Vries, 1958 [8] ); L w is the latent heat of water vaporization (J m −3 ); λ(θ) is the soil thermal conductivity (W m -1 K -1 ); and q L and q v are the liquid water flux and vapor flux, respectively (m s -1 ).
Neglecting hysteresis between drying and wetting curves, the van Genuchten model (van Genuchten, 1980 [9] ) is adopted to simulate the Soil Water Characteristic Curve (SWCC), as follows:

( )
where Θ is the effective saturation (dimensionless) and Θ=(θ L -θ r )/(θ s -θ r ), θ s and θ r are the saturated and residual water content (m 3 m -3 ), respectively; α (m -1 ), n (dimensionless) and m (=1-1/n) are empirical shape parameters.The isothermal unsaturated hydraulic conductivity of liquid water in both frozen and unfrozen soil is modeled following Mualem (1976) [10] : where K s is the saturated hydraulic conductivity (m s -1 ); l is a fitting parameter and was set to 0.5 by Mualem (1976) [10] .The definitions of thermal hydraulic conductivity of liquid water K LT , and hydraulic conductivities of vapor due to temperature K vT , and matric head K vh refer to Saito et al. ( 2006) [11] : respectively.In Eq. ( 5), G wT is the gain factor which amends the temperature dependence of the SWRC; γ is the surface tension of soil water (J m -2 ); and γ 0 is the surface tension at 25 ºC (= 71.89 g s −2 ).In Eqs. ( 6) and ( 7); M is the molecular weight of water (= 0.018 kg mol - 1 ); g is the gravitational acceleration (m/s 2 ); R is the universal gas constant (= 8.341 J mol -1 K -1 ); H r is the relative humidity (dimensionless); η is an enhancement factor (dimensionless); ρ vs is the saturated vapor density (kg m -3 ); and D is the vapor diffusivity in soil (m 2 s -1 ).
The parameter H r relates vapor density to matric head and temperature and is expressed is as follows: exp Definitions of other parameters used in this model can be found in Teng et al. (2015) [7] .

Numerical analysis for canopy effect
Equations ( 1) and ( 2), subjected to appropriate initial and boundary conditions, are solved numerically using the finite element method for spatial discretization and the finite difference method for temporal discretization.The solution procedure is facilitated by the COMSOL Multiphysics package (5.0).In order to reveal the mechanism of canopy effect, two cases are designed for computation as presented in Table 1.A one-dimensional soil column is chosen for the analysis, with a height of either 5 m or 20 m.The water table is located at the bottom of the soil column, while the top boundary either has an evaporating flux (10 cm/d) or is completely sealed (zero flux).The initial water content is set to 10% and is uniform over the depth.The initial temperature along the soil profile is 15 °C.The top and bottom boundary temperatures are constantly maintained at 0°C, and 15 o C, respectively, for a period of 90 days.
Considering the engineering background of canopy effect as stated in Li et al. (2014) [1] , this study chose the silty soil in Lanzhou area for analysis.The Soil Water Characteristic Curve of the silty soil is given in Figure 2, which was simulated by the van Genuchten model.Figure 3 presents the simulated unsaturated permeability by the Mualem model.The computed results of case (1) are shown in Figure 4.
The figures numbered (a) and (b) denote the groundwater table located at 5 m and 20 m, respectively.The impervious cover affected the moisture content to a depth of 50 cm when the groundwater table is located at 5 m, and a depth of 100 m when the groundwater table is located at 20 m.The maximum difference between the covered and uncovered moisture profiles is about 4% of volumetric water content, which takes place at the ground surface.It is also interesting to note that the moisture content underneath the cover does not increase above the initial value of 10%, which is above the residual moisture content.The capillary rises in Figure 4(a) and 4(b) are around 4.5 m.The impervious cover seems to stop the evaporation at the ground surface, but fails to accumulate significant moisture.The canopy effect reported in Li et al. (2014) [1] is not apparent in this case.Similar results shown in Figure 4 have also been obtained using commercial software such as Hydrus-1D [12] .Test runs using the COMSOL model and Hydrus-1D fail to reproduce a significant amount of moisture accumulation under impervious covers, which implies that existing models based on the PDV theory is incapable of generating the canopy effect of type A 2 .A closer inspection of the PDV theory reveals its shortcomings.According to the PDV theory, liquid water flow and vapour flow is driven either by a matric suction gradient or by a temperature gradient.The vapour density is related to the matric suction and temperature through Eq. ( 8), while the matric suction is related to the moisture content via the SWCC (Eq.( 3)).If evaporation is stopped by an impervious cover, or vapour condenses into liquid water, the liquid water content in the soil increases and the matric suction consequently decreases, which in turn increases the vapour density and reduces the moisture flow.Therefore, existing models based on the PDV theory predict an equilibrium state where moisture flow driven by temperature gradients in unsaturated soils is balanced by that driven by the suction gradient.Such a case is very similar to a cup of hot water covered by a cold ceramic lid.Some vapour condensation will take place at the cold ceramic lid, but not much water will disappear from the cup, as the evaporation will stop when there is sufficient liquid water under the lid.However, the situation will change if the ceramic cover is replaced by a cover that absorbs water, for example, a cover made of dry wood.The water in the cup will transfer to the absorptive cover until it is fully saturated.The key question is then: what is the absorption mechanism for a moisture profile like A 2 shown in Figure 1? Or alternatively, under what mechanism can the moisture content continue to increase without increasing vapour density or decreasing matric suction?

A new mechanism for canopy effect
We first note that the two cases reported in the literature with significant moisture accumulation underneath impervious covers are both related to frozen soils (Li et al. (2014) [1] ; Eigenbrod and Kennepohl (1996) [2] ).
Another inspirational example is that ice tends to grow continuously inside an old type freezer (the non-frost-free type) as long as warm and humid air enters the freezer.Realising these facts, we propose here a new mechanism of canopy effect as follows.
When the temperature in a soil drops below the frost point, vapour under an impervious cover will change into ice through desublimation.The formation of ice will further reduce air humility and increase the matric suction, which in turn accelerates the vapour transfer in the soil.In other words, the vapour-ice desublimation acts like a sink term in the moisture transfer equation, and releases energy (latent heat) in the heat transfer equation.This process could lead to significant accumulation of total water content (including ice) under the cover, and possibly a moisture profile of type A 2 .The terminology 'canopy effect' is justified because of the absorption effect of the vapour-ice desublimation process.
A specific laboratory experimental program is designed to validate the above mechanism.The test results strongly support the freezing induced canopy effect.However, due to the page limitations, the experimental program and the test results will not be presented here.Instead, this paper focuses on the numerical simulation of the phenomenon.

Theory for vapour transfer in unsaturated freezing soil
Unsaturated soil mechanics so far mainly applies to soils with super-freezing temperatures.Thus, a model for coupled liquid water, vapour, and heat movements in unsaturated soils subjected to both sub-and superfreezing temperatures is not available in the literature.Such a model is proposed here, to reproduce the canopy effect observed in field.
In a sub-freezing condition, water exists in soil pores in three forms: liquid water, water vapour and ice.The endothermic and exothermic processes are accompanied with phase changes.Here, we extend the governing equations for coupled heat and mass transfer, i.e.Eq. ( 1) and ( 2), to include phase change.The mass conservation is expressed as, where θ i is pore ice content (m 3 m -3 ); Kʹ Lh is the permeability in unsaturated freezing soil (m/s) and differs from the permeability of unsaturated soil K Lh .The expression of Kʹ Lh refers to Taylor and Luthin (1978) [13] : The adjusting parameter 10 -10θi denotes the obstruction of pore ice to vapour flow.When the soil temperature is over 0°C, θ i becomes zero, and Kʹ Lh reverts into K Lh as in Eq.( 3).
The energy balance in an unsaturated freezing soil is described as follows, The four terms on the right-hand side of Eq. ( 11) are the soil heat flow by conduction, the convection of sensible heat due to flowing water, the transfer of latent heat by vapor diffusion, and the transfer of sensible heat by vapor diffusion, respectively.In Eq. (11), λʹ(θ) is the heat conductivity of unsaturated freezing soil (W m -1 K -1 ), and is obtained by replacing θ L of λʹ(θ) in Eq. ( 2) with θ L +Fθ i (Hansson et al., 2004) [14] .
As for the movement of moisture and heat in the frozen zone, an additional constitutive relation is required to solve the governing equations.Here, the generalized Clapeyron equation is used to establish the relationship between the pressure of unfrozen liquid water and that of pore ice when these two phase are in equilibrium (Sheng et al., 2014) [15] : where T 0 is the freezing point of water (K); and T in this equation is in unit of ºC.

Numerical reproduction of canopy effect
The proposed theory was numerically solved by using COMSOL Multiphysics (5.0), with vapour transfer and ice formation taken into account.In attempt to numerically reproduce the canopy effect using the proposed theory, we design a numerical example as follows.A soil column of 20 m in height is subjected to a temperature of -5 °C at its top and a temperature of 15 °C at its bottom.The top boundary of the soil is sealed for any water flux.The initial water content is 7%, and the computation period is 90 days.
Figure 5 shows the computed profiles of total water content (θ L +θ v +θ i ), temperature, liquid water flux and vapour flux along the soil depth.Figure 5(a) depicts that the water content near the lower end is significantly raised from the initial value, with a capillary acting between 16-20 m at 90 days.The water content between 2.5 m and 16 m stays close to the initial level.For the soil within the top 2.5 m, the water content shows a significant increase, especially for the soil at top 1.0 m where a saturated water content of 42% is reached.Figure 6(b) gives the depths of the freezing front (0°C line), which is located at about 0.2 m, 0.5 m, 0.8 m respectively for 5 days, 30 days and 90 days.These positions coincide with the depths of maximum water content in Figure 5(a).Moreover, Figure 5(b) shows that the temperature profile below 8 m remains constant during the freezing period, which agrees with the real condition in Lanzhou area.In Figure 5(c), the positive value represents the upward flow of liquid water, and negative downward.The positive value of water flux between 16-20 m means that the liquid water in unsaturated soil migrates upwards under the capillary force, and the water content profile reaches an equilibrium state after a certain period.For the soil between 2.5-16 m, little water flux occurs there.The negative value of water flux in the top 2.5 m soil indicates a small downward liquid flow.The reason might be that the unfrozen water content near the surface is higher than that at a deeper position.It also implies that the increase of water content within the top region is not caused by the liquid water flow.The embedded chart in Figure 5(c) shows that the liquid water flux is close to zero within a top layer of soil (0-0.2 m at 5 days, 0-0.5 m at 30 days and 0-0.8 m at 90 days), which is due to the phase change of water into ice.Figure 5(d) presents the computed results of vapour flux, and shows a clear upward flux near the surface zone.It can be concluded that water accumulation beneath the impervious surface is caused by vapour transfer.The vapor flux reaches a peak value at a certain depth, which is likely due to the vapour condensation and ice formation that occupy the soil pores and impede vapor transfer.
It is interesting to note that moisture would not accumulate in an unsaturated freezing soil exposed to atmosphere.In this case, evaporation and sublimation at the ground surface prevents ice accumulation in the soil, a process analogous to the mechanism in a modern frostfree freezer.Inside such a freezer, air is circulated by means of one or more fans, and air circulation helps sublimate any ice or frost that may form in the freezer.Another example is that wet clothes hanging outside in cold winter will eventually dry up, even though some ice may initially form in the clothes.Numerical simulation of an unsaturated freezing soil with an open top surface is difficult, as the boundary condition of evaporation and sublimation is not easy to specify.However, field observation shows that no ice forms in very dry soils that are open to the atmosphere.
The canopy effect is reproduced in the numerical results in Figure 5, which shows that the canopy effect is an outcome of the interaction between vapour transfer in the unsaturated soil, the ice formation in the freezing soil and the impervious cover.The temperature gradient drives up some moisture flow towards to the top soil, but is not sufficient to saturate the soil, particularly when the groundwater table is low.Formation of ice in the soil further increases the matric suction and decreases the vapour density, which accelerates the upward vapour flow, leading to significant total moisture accumulation.
The model presented in this paper is one-dimensional only.In reality, vapour transfer in soil can take place in multiple directions.Humid air may enter the vapour-ice desublimation zone from horizontal directions, further enhancing the canopy effect.

Conclusions
The role of vapour flow is usually ignored when analysing water migration in unsaturated soils.Indeed, vapour flow is seldom considered in geotechnical design.However, there are circumstance where vapour transfer can lead to engineering problems and needs to be considered.One of such circumstances is the moisture accumulation caused by vapour transfer under impervious covers, and is referred to as the canopy effect.It was reported that significant moisture could accumulate in unsaturated freezing soils, which led to frost damages to airport runways and road pavements.
Existing models based on the Philip and de Vries theory of coupled heat and mass transfer in unsaturated soils fall short of capturing significant canopy effect.
A new approach of modelling the coupled movements of moisture and heat in unsaturated freezing soils is presented in this paper, in which phase changes of evaporation, condensation and de-sublimation are taken into account.The computed results demonstrate that soil freezing fed by vapour transfer can produce a water content close to full saturation for the soil under an impervious cover.Systematic experimental verification of the proposed numerical model is warranted.

Figure 2 .Figure 3 .
Figure 2. SWCC of silty soil at Lanzhou area, with symbols representing measured data, solid line the van Genuchten model.

Figure 4 .
Figure 4. Moisture content profiles at 90 day, with an initial moisture content of 10%, and a constant water table at (a) 5m, and (b) 20 m deep, respectively.