Influence of the soil-atmosphere exchange on the hydric profile induced in soil-structure system

Soil-atmosphere exchange leads to a moisture change in the soil. This can cause major damage to engineering structures due to the soil expansion and shrinkage. The soil-atmosphere exchange is related to several parameters, in particular the soil characteristics and climate conditions. The presence of an engineering structure causes a variation of the hydraulic profile in the soil, which can lead to heterogeneous soil movement and consequently to structural damage. This paper presents a coupled numerical model based on the consideration of both water flow in unsaturated soils and soilatmosphere exchange. After the validation of the model, the paper presents its use for the analysis of the influence of the presence of structures on moisture change induced under climatic conditions recorded in a semi-arid region. Analysis shows that the presence of the structure leads to important change in the moisture distribution, in particular in the vicinity of the structure.


Introduction
The interaction between the atmosphere and soil could induce a significant change in water saturation in soils with a risk of damage to engineering structures.This damage, especially for lightweight structures, is often observed after cycles of swelling and drying (Blight, 1997).During a drought period, the soil water content in the vadose zone is reduced because of evaporation.Settlement can occur as a result of suction increase and soil shrinkage.It should be noted that two principal processes govern the exchange of water between the soil surface and the atmosphere.Water enters the soil surface as liquid through the process of infiltration.Alternately, water exfiltrates from the soil surface as vapor through the process of evaporation.The process of infiltration depends primarily on soil properties such as hydraulic conductivity and is reasonably well understood.The evaluation of the evaporative fluxes from a soil surface is more diffi-cult, since the rate of evaporation depends on both soil properties and climatic conditions (Wilson et al., 1994).The coupled numerical models of the soil moisture and heat balance in porous media may be used to describe the soil-atmosphere exchange.Combination of these models with the equation of the lower atmosphere boundary layer enables the estimation of evaporation and energy balance at the soil surface.Attempts have been made in recent years to solve the coupled equation of heat and moisture transfer using the FEM technique.Wilson and Frdlund (2000) provided the results of SoilCover program, which takes into account changes in climate conditions.Yakirevich (1997) presented a mathematical model using a finite element code with an implicit scheme for calculating the evaporation of saline soils considering the migration of salt.Mangeney et al. (2003) developed a code SiSPAT, which takes into account several soil layers.This code has been used mainly in agricultural applications.Yanful and Mousavi (2003) used the SoilCover code to study the influence of the presence of several soil layers in conserving the moisture in the soil.Yanful et al. (2003) estimated the evaporation on the soil surface and compared the results with those of Wilson (1990).Gitirana et al. (2005Gitirana et al. ( , 2006) ) employed a two-dimensional model to study the influence of soil-atmosphere exchange on slope stability.Vu et al. (2007) developed a 2-D model to study the influence of climatic condition and soil properties on evaporation.This paper presents a numerical analysis of the influence of the presence of a light structure such as a house building on the hydric profile induced by the soil-atmosphere exchange.This analysis allows the determination of the suction induced by drought around the structure, and then the evaluation the drought-induced soil movement which could cause structure damage.Analysis is carried out using a coupled numerical model based on the consideration of both water and heat flow in unsaturated soils and the soil-atmosphere exchange.Since this analysis concerns light structure, the variation in the soil properties due to the structure construction is not considered in the numerical model.Analysis shows that the soil-atmosphere exchange induces a significant variation in the suction around the structure lateral boundaries, in particular for clayey soils.This variation could lead to heterogeneous soil movement, and consequently to structural damage.

Mathematical formulation
The governing partial differential equations for coupled moisture and heat flow in soils are derived using the principle of conservation of mass and energy.In general, the principle of conservation may be stated as ∂u ∂t = −∇q where u is the amount of mass or heat per unit volume, q is the local average flux and t is the time.The symbol ∇ denotes the Laplace operator.As a consequence of evapotranspiration, a transient water flow takes place through the soil from the water table to the soil surface.Heat flow simultaneously occurs because of the temperature gradient.Water flow, heat flow and volume change are fully coupled.For simplicity in multiphase modeling, Wilson et al. (1994) considered the soil to be rigid.Under constant net stress conditions, Wilson et al. (1994) considered the following expression relating the change in water content to the suction change; The mass transfer equation can be derived directly from the Richards equation for transient flow in unsaturated soils with adaptation for vapor flow added by Wilson (1990).The driving force for Darcian liquid flow is the total head but the equations have been re-written as functions of pressure so that the hydraulic heat and mass equations can be fully integrated.
The generalization of the 1-D Wilson's equation (1990) for the transfer of liquid water and water vapor in an unsaturated 2-D soil can be written as follows: where ρ is the density of water, P v is the vapor pressure of the soil moisture, H is the hydraulic head in the water phase (m) i.e.,u w ρg + y , u w is the pore-water pressure, y is the spatial position, g is the acceleration due to gravity, K x and K y are the hydraulic conductivity of soil in the x and y directions, respectively.It should be noted that D v is the vapor diffusion coefficient described by Wilson (1990) [ , where D vap is the coefficient of water vapor diffusion in free air.According to Kimball et al. (1976) (Wilson, 1990).(Gitirana et al., 2005).
stands for the temperature; m a is the molecular mass of water vapor (m a = 18.016 kg/kmol), R is the universal gas constant (R=8.314j/mol.K), β=(1 − S r )n is the cross sectional area of the soil available for the vapor flow; S r is the degree of saturation and n is the porosity.α=β 2/3 is the tortuosity factor and m w is the slope of the volumetric water content function.
By coupling the equation of the flow of steam given by Fick's equation and the heat transfer given by Philip and de Vries (1957), we obtain the following equation for the heat transfer in the soil: Where K tx and K ty are the thermal conductivity of the soil in the x and y directions, respectively; L v is the latent heat of vaporization of water, which is a function of temperature as proposed by Bertin et al. (1981) (L v = 4186×(607-0.7×T )); λ t is the volumetric specific heat.Equations ( 2 numerical procedures.This drawback can be overcome using the equation (Edlefsen and Anderson, 1943) : where P vs is the saturated vapor pressure of pure free water, S is the suction and h r is the relative humidity.Joshi (1993) reformulated the above equation as follows: (5) d 1 and d 2 are given by the following expressions: In order to simulate natural systems, the moisture and heat flow equations have to be coupled to the atmosphere through evaporation, precipitation, and heat flux terms.This requires evaluating the boundary conditions at the soil-atmosphere interface.Wilson (1990) modified the equation of Penman (1948) for the determination of the evaporation (E vap ) at the surface of unsaturated soils as follows: is the slope of the saturation vapor pressure versus temperature curve at the mean temperature of the air, R n is the net radiation at the soil surface, L v is the latent heat to vapor, η is the psychrometric constant, U a is the wind speed, P a is the vapor pressure in the air above the evaporation surface; φ a and φ s designate the relative humidity in the air and at the soil surface, respectively.
The surface temperature may be estimated using the Wilson's equation (1990): T s is the temperature at the soil surface and T a is the temperature of the air above the soil surface.
The following equation is used for permeability (Fredlund et al., 1994;Leong and Rahardjo, 1997): K s is the saturated permeability, a is the air-entry value of the soil and n,m and p are materials constants, e is equal to 2.718.
For the water retention curve, the Van Genuchten equation ( 1980) is used: where θ is the volumetric water content, θ r is the residual volumetric water contents, θ s is the saturated volumetric water contents, and a s ,n s and m s are material constants.
For the thermal conductivity, K t , the equation given by Van de Griend and O'Neill (1986) is used: s is a texture-dependent coefficient; C sec stands for the volumetric specific heat capacity of the soil solids (4.18×10 6 j/m 3 K); The coefficient of the specific heat is provided by de Vries (1963): Where λ h is the volumetric specific heat, C sec is the volumetric specific heat capacity of the soil solids (2.24×10 6 j/m 3• C), C w is the volumetric specific heat capacity of the liquid water phase (4.15×10 6 j/m 3• C), C a is the volumetric specific heat capacity of the air phase, which can be assumed to be negligible.θ s , θ w , θ a volumetric contents of solid, water and air, respectively.

Numerical implementation
The soil atmosphere-exchange model presented in the previous sections has been implemented in the finite element code ESNA, which has been developed at the Laboratoire de Mécanique de Lille.After substituting Eq. ( 4 and (3), the system of the ordinary equations can be obtained on the current configuration as follows: B is the gradient of the interpolation matrix, M is the mass matrix and S and T are the nodal values of suction and temperature, respectively.By applying a time integration scheme, one can obtain the algebraic counterparts of Eqs. ( 14) and ( 15).Assuming that the values of suction, temperature and their derivatives { S t , Ṡt ,T t , Ṫ t } are known at time t, the integration consists of updating {S t+ t , Ṡt+ t ,T t+ t , Ṫt+ t } at the next step t+ t, using to the Crank-Nicholson implicit scheme as follows: K w , K wt , K tw and K t represent the global matrices containing the coefficients which appear in Eqs. ( 14) and (15).F w and F t denote the heat and flux vectors, respectively.It should be noted that for the developed soil atmosphereexchange model, two different methods are implemented to determine the soil parameters.In the first method, the material parameters K,K t , λ t and m w are interpolated for each time step using available experimental data.In the second method, the material parameters are calculated from Eq. ( 10), ( 11), ( 12) and ( 13  values of temperature and suction in order to evaluate the evaporation and heat flux at the ground surface as well as the soil parameters required for the determination of the current temperature and suction values.

Verification of the numerical model
The FEM code was validated using the laboratory test of Wilson (1990) and the numerical simulation of this test conducted by Gitirana (2005).Wilson (1990) performed a column-drying test in which the change in temperature and volumetric water content were monitored.The column was filled with saturated Beaver Creek sand, and was allowed to drain at the bottom so as to reach a hydrostatic condition.The column was then sealed at the base and at the circumference.The exposed surface of the sand was allowed to evaporate under controlled laboratory conditions with a constant temperature of 38 • C and a relative humidity of approximately 10 % for a period of 35 days.The temperature profile was monitored using 12 thermocouples installed in the sand at different depths.The water content profile was measured on soil samples extracted from 12 ports at different levels.The drying column experiment and its initial and boundary conditions are shown in Fig. 1.
In this test, water and heat flow at the bottom and lateral boundaries are zero.The upper boundary is in contact with the atmosphere.Figure 2a presents the hydraulic soil properties of the soil.Figure 2b presents the thermal property function of the Beaver Creek sand.These functions were obtained using the formulations presented by de Vries (1963).According to Eq. ( 13), the thermal conductivity and the volumetric specific heat are functions of the characteristics of the individual phases and functions of the amount of water stored in the soil as well.The effect of the amount of water stored in the soil pores is shown by the decrease in both thermal conductivity and volumetric specific heat, as soil suction increases.Figure 3a and b shows a comparison of the numerical results with both the experimental results of Wilson (1990) and the numerical simulation of Gitirana (2005).As it can be seen, the numerical results are in good agreement with those of Wilson (1990) and Gitirana (2005).It can be concluded that the proposed numerical model is capable of reasonably predicting both water content and temperature profiles on the basis of measured atmospheric boundary conditions and soil properties.

Case study
The numerical model is used to study the influence of the presence of an engineering structure on the soil hydraulic profile, which is induced by the soil-atmosphere exchange.
Figure 4 illustrates a general view and dimensions of the case study.The study concerns a light structure supported by a raft foundation.In the numerical model, only the raft foundation is considered, because the substructure does not affect the soil-atmosphere exchange.In addition, since the construction of light structures does not significantly disturb the soil material, the change in the soil properties due to the structure construction is not considered in the numerical model.The depth of the soil layer is considered as 3 m, while the width of the model is found from a sensibility analysis to be 30 m after which the structure has a negligible effect on the suction profile.The base and sides of the soil are considered to be impermeable.Exchange conditions are applied on the upper surface.Figure 5a shows the water retention curve that is given by Eq. ( 11) whose parameters are given in Table 1.The permeability, according to Eq. ( 10) is shown in Fig. 5b, whose parameters are given in Table 2. Figure 5c and  d shows the variation of the thermal conductivity and specific heat with water content, respectively.The thermal conductivity K t is given by the formula of Van de Griend (1986) which is described in Eq. ( 12) whereas the specific heat, λ t , is obtained from Eq. ( 13) (de Vries, 1963).The climatic conditions used in the analysis were recorded in June 2005 in a semi-arid area in the South of Syria (Fig. 6).

Analysis in the absence of the engineering structure
The case study is first analyzed in the absence of any engineering structure.Figure 7a shows the variation of the  evaporation rate.It can be seen that the general trend of evaporation rate decreases with time.The low value of evaporation at t = 15 days is due to the relatively high precipitation (Fig. 7a).The suction variation at the top of the soil layer is shown in Fig. 7b.It can be seen that after five days the suction slightly increases to 10 kPa.Then a sharp increase in suction is observed.After 10 days, it attains 3 MPa.The sharp increase in suction results from the decrease in the permeability over the same period.After the sharp increase, the suction augments steadily except during precipitations, which cause a suction decrease.Figure 7c shows the computed permeability at three different locations.It can be seen that after a slight decrease over six days, there is a large decrease in the permeability over the following four days.This decrease leads to a reduced water flow towards the soil surface, which explains the large increase in the suction at the soil surface over the same period.Figure 7d shows the suction profile in the soil.It can be observed that the suction increases, in particular in the vicinity of the soil surface.After 30 days, the 28 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8  suction augments up to 1 MPa at depth z = 1.1 m and only to 100 kPa at z = 1.75 m; thereafter the suction remains approximately constant in the range 40-70 kPa.

Influence of the presence of an engineering structure
The engineering structure is considered as an impermeable border (i.e., soil-atmosphere exchange is not permitted).The width of the structure is assumed to be 5 m (Fig. 8). Figure 9a shows the computed suction profiles at different points in the soil mass.It can be seen that the variation of suction in point A, which is situated about 25 m from the structure, is quite similar to that calculated in the previous case without the structure.This shows that the model boundary is located sufficiently far from the foundation.At point B, which is situated in the vicinity of the structure, a significant increase in suction occurs after 8 days.At point C, which is situated underneath the centre of the structure, a very small increase in the suction is observed.Figure 9b shows the distribution the suction at the soil surface at t = 5, 15 and 30 days.It can be noted that suction distribution undergoes a sharp variation at the junction between the soil and the structure.After 30 days, the suction changes from approximately 100 kPa underneath the structure to about 100 MPa along the free surface.

Protection of engineering structures against soilatmosphere exchange
Protection of the structure against the drought can be ensured by the construction of an impervious area in the vicinity of the structure (Fig. 10).To study the effectiveness of this technique, simulations were performed with three values of the length of the protection area (LC = 1, 1.5 and 2 m). Figure 11a shows the influence of the protection area on the suction variation under the foundation (point C) and at board (point B).It can be noted that the presence of the protection area reduces the development of the suction in the vicinity of the foundation significantly.Figure 11b shows the variation of suction at the soil surface.This confirms the observation in Fig. 11a: The construction of a protection reduces the gradient of suction under the foundation and therefore the differential settlement induced by the drought.The width of the protection area affects the extent of the zone of quasiuniform distribution of the suction under the foundation.

Influence of the type of soil (parametric analysis)
The type of soil plays an important role in the phenomenon of transfer.In this section, we analyze its influence on the suction induced by the soil-atmosphere exchange in the vicinity of the structure.Calculations were performed with two soil types (clay and silt) whose parameters are given in Tables 3  and 4. For the same level of suction, the permeability of silt is higher than that of clay (Fig. 12a), while the water content of clay is higher than that of silt (Fig. 12b).Figure 13a and b shows the influence of the soil type on the variation of the suction and permeability at point A. For the  silty soil, the suction increases during the observation period up to 100 kPa, whereas in the clayey soil, there is a sudden increase in suction at the 6th day reaching about 100 MPa.
Figure 13b shows that the permeability of the silt undergoes a smaller variation than that of the clay.
Figure 13c shows the influence of the soil type on the variation of the suction in a vertical section of the soil mass.For the silty soil, a uniform distribution of suction with depth is observed.For the clayey soil, we observe a sharp variation of the suction in the surface layer.This variation is due to the drastic reduction of clay permeability, which leads to a concentration of the mass transfer near the soil surface.
Figure 13d shows the influence of the soil type on the variation of suction at the soil surface.With the clayey soil, we observe a significant difference between the values of suction under the foundation and the free surface.For the silty soil, the suction under the foundation increases and remains close to that at the free surface.This can be explained by the presence of a lateral water transfer in the silt.

Conclusions
This paper presented an analysis of the influence of the presence of a light structure on the soil hydraulic profile induced by the soil-atmosphere exchange.The analysis was conducted using a coupled finite element modeling that integrates both the flow in unsaturated soil and the soil-atmosphere exchange.This model was used for studying the influence of the soil-atmosphere exchange in a semi-arid area on the soil hydric profile.Analyses show that the soil-atmosphere exchange induces uniform distribution of suction with depth in silty soils, but a high concentration of suction in the surface layer in clayey soils.This concentration is due to the drastic reduction of clay permeability in the surface layer.Numerical results show that the presence of a structure affects the distribution in the vicinity of the structure of the suction induced by climatic variation significantly.In clayey soils, the presence of the structure induces a sharp variation of suction in an area located around the lateral boundary of the structure.This sharp variation is expected to induce a differential settlement with eventual structure damage.The construction of an impervious area around the structure reduces the lateral variation of the moisture under the structure significantly.It could effectively reduce the differential settlement and consequently reduce the risk of structure damage due to the soil-atmosphere exchange.
Edited by: F. Castelli Reviewed by: two anonymous referees

Fig. 4 .Figure 5b .
Fig. 4. Configuration used to analyse the influence of the presence of the structure on the soil-atmposphere exchange.

Fig. 6 .
Figure 6.Atmospheric condition used in the case study (Recorded in June 2005 in the south of Syria) Fig. 6.Atmospheric condition used in the case study (Recorded in June 2005 in the South of Syria).

Figure 7b .Figure 7b .
Figure 7a.Variation of the evaporation at the soil surface (Case study)

Fig. 7c .
Figure 7c.Variation of the soil permeability

Figure 7d .
Figure 7c.Variation of the soil permeability

Figure 9a .Figure 9b .
Figure 9a.Influence of the presence of the structure on the variation of soil suction

Fig. 9a .Figure 9b .
Fig. 9a.Influence of the presence of the structure on the variation of soil suction.

26Figure 10 .Fig. 10 .
Figure 10.Numerical modelling of the protection area around the structure Fig. 10.Numerical modelling of the protection area around the structure.

Figure 11a .Figure 11b .Fig. 11a .
Figure 11a.Influence of the protection area on the varition of the soil suction

Figure 11a .Figure 11b .
Figure 11a.Influence of the protection area on the varition of the soil suction

Figure 12a .Figure 12b .
Figure 12a.Hydraulic conductivity of the clay and silt considred in the parametric study

Figure 12a .Figure 12b .
Figure 12a.Hydraulic conductivity of the clay and silt considred in the parametric study

Figure 13a .Figure 13b .Fig. 13a .
Figure 13a.Influnce the soil type on the variation of the suction in point A

Table 3 .
Influence of the soil type; retention curve parameters for clay and silt.

Table 4 .
Influence of the soil type; permeability parameters for clay and silt.