Physical and numerical modelling of soil-atmosphere-structure interaction

. Extreme, extended wet and dry seasons increase the adverse effects that soil wetting and drying cycles have on the response of shallow geotechnical structures. In expansive soils, volumetric changes due to water content variations may result in the incompatibility of deformations at the soil-structure interface. This study proposes a physical approach and a numerical model to address the soil-atmosphere-structure interactions during soil saturation and desiccation. Experimental desiccation tests were performed on relatively thin, compacted kaolin clay samples that represent the soil-atmosphere boundary. A climatic chamber was used to impose atmospheric conditions of air relative humidity, temperature, wind velocity, and irradiance on the soil surface. Empirical mathematical expressions were obtained to estimate soil desiccation rates as a function of basic atmospheric parameters and soil properties. The experimental desiccation approach was complemented with a coupled thermo-hydro-mechanical (THM) numerical model for unsaturated soils. The coupled THM model calculates water and thermal fluxes, soil volumetric changes, vertical stresses, and the compatibility of soil-structure movements during swelling and shrinking. An example of the capabilities of the numerical model is presented for a full-scale geotechnical problem.


Introduction
The effects of climate change on the analysis, design, construction, and rehabilitation of geotechnical projects have become of great importance as more frequent extreme, extended wet and dry seasons occur in different places of the world. The performance of structures supported near the soil surface depends on soil and environment interaction. For instance, shallow geotechnical structures constructed on expansive soils are subjected to soil volumetric changes due to wetting and drying cycles (i.e., swelling and shrinkage) that may result in excessive deformations and instability.
Understanding the interaction mechanisms between soil atmosphere and soil structure is critical in evaluating the potential effects of soil movement on these shallow structures. Two key soil-atmosphere interaction processes govern surficial expansive soils' behaviour and mechanical response. Water infiltration produces heave movements and upward pressures beneath foundations. On the contrary, water evaporation produces desiccation, shrinkage, and cracking. Both upward and downward movements affect the stability of shallow structures.
The water infiltration process in unsaturated porous media has been widely studied and calculated with relatively good levels of accuracy [1,2,3]. However, soil evaporation for geotechnical engineering purposes is a relatively more complex phenomenon due to the challenges in coupling the atmospheric conditions to the unsaturated soil properties [4,5,6]. The complexity of the analysis for both the wetting and drying scenarios increases when the structural response is assessed.
Several atmospheric conditions, soil properties, and structural considerations must be considered in the geotechnical analysis of soil-atmosphere-structure interaction problems. The main atmospheric conditions at the soil-atmosphere interface include precipitation, air relative humidity, air temperature, wind velocity, solar radiation, and atmospheric pressure. The main soil properties involved in water and heat fluxes and in volume changes include hydraulic and thermal conductivity, porosity, water retention capacity, and stiffness. From the structural standpoint, properties of the foundation system such as the geometry, structural stiffness, and loading conditions determine the concentration of stresses, structural capacity, and allowable deflections. The interactions among the above-mentioned properties and atmospheric conditions need to be coupled in a thermo-hydro-mechanical (THM) model that simultaneously calculates the thermal, hydraulic, and mechanical response of the unsaturated porous media and the structure.
This study is divided into two components. The first part focuses on soil desiccation rates. It describes the experimental approach used to estimate evaporation rates for a wide range of atmospheric conditions. The results are used in the numerical model. The second part summarises the numerical approach used to model soilatmosphere-structure interactions. A THM model using the Explicit Finite Difference Method (E-FDM) is proposed. This section includes the analysis of a singlestory building supported on a plastic soil.

Soil-atmosphere interaction
The study of soil-atmosphere interactions addresses the physical and phenomenological relationships that exist between the soil's hydraulic and thermal properties and the atmospheric conditions involved in the water balance (i.e., infiltration, runoff, evaporation, and evapotranspiration). In a simplified manner, water evaporation from soils occurs when the atmospheric conditions above the soil surface are drier than the soil surface. If the atmospheric conditions are inverse, the soil gains water. Hydraulic and thermal transfer interactions have place at the soil-atmosphere interface. If the transfer mechanisms (i.e., water, vapour, and thermal exchanges) at this interface are well understood, the soil surface becomes a suitable boundary condition for soil discretisation and modelling.
The experimental tests in this study focused on evaporation from the soil-atmosphere interface. Infiltration was not experimentally determined but was physically accounted for in the modelling section. Evapotranspiration and runoff are not part of this study.

Soil-atmosphere-structure interaction
Geotechnical engineering problems deal with the soil response to loading and climate. In near-surface, lightlyloaded geotechnical structures such as shallow foundations, the effects of climate may be as important as the soil structural capacity. Soil shrinking and swelling in expansive soils may result in total or differential movements beyond the serviceability limits of the structure. In fact, differential movements between the soil and the foundation may develop void zones beneath the foundation during shrinkage and uplifting during swelling. These movements increase the concentration of stresses on the structural and nonstructural members. A coupled THM model for analysing soilatmosphere-structure interaction problems needs to address the thermal and hydraulic exchanges at the soilatmosphere boundary and the compatibility of stressstrains between the soil and the structure.

Experimental evaporation testing
In soils, evaporation starts at the soil-atmosphere interface and extends to variable depths below the ground surface. This is known as a drying or evaporation front, and its depth and shape depend on the properties of the soil (i.e., water retention capacity), groundwater level, and intensity of the atmospheric conditions at the soil-atmosphere boundary (i.e., the evaporative energy produced by the atmospheric conditions). Given the importance of the soil-atmosphere interface on the water, vapour, and thermal fluxes, this study evaluated the evaporation process in relatively thin soil samples representing the interface boundary. A brief description of the climatic chamber and soil used in this study is provided as follows.

Equipment
Evaporation tests were performed using the climatic chamber designed and constructed by the Universidad de Los Andes. The climatic chamber is a robust equipment that simulates dry atmospheric conditions on the surface of relatively thin soil samples. Soil specimens are accommodated on glass containers of various heights. Insulation materials (i.e., expanded polystyrene foam) are placed along the perimeter of the glass container and beneath the glass to decrease heat fluxes from the sides and bottom of the soil. The chamber operates with a centrifugal fan that induces the airflow above the soil sample, a heating system to increase the air temperature to the target temperature, a set of infrared lamps to simulate irradiance on the soil surface, and a cooling system to dehumidify the air. The walls of the chamber are made of plexiglass. The schematic plan view of the climatic chamber is shown in Fig. 1. Various types of sensors are installed on a moving electronic module above the soil surface to measure the atmospheric conditions during testing. Wind velocities are measured with a couple of air velocity sensors. Relative humidity is measured with a set of relative humidity sensors. Thermocouples are used to measure air and soil temperatures. The soil surface temperature is measured with an infrared sensor. The evaporated mass is measured with an electronic scale. Data is collected with an acquisition system. The schematic 3D view of the climatic chamber is shown in Fig. 2. Details on the design and operation of the climatic chamber may be found in [7].

Materials
Desiccation tests were performed on a plastic, compacted kaolin clay compacted to its optimum moisture content and maximum dry density as per the standard compactive effort. A compacted material was preferred over a natural soil to decrease the variability of the soil properties. Soil sample thicknesses varied between 3 and 20 mm. The soil surface area for all the specimens was 240 x 240 mm 2 . The main index properties, optimum moisture content (OMC), and maximum dry unit weight of the soil are summarized in Table 1. The main hydraulic, thermal, and mechanical properties of the soil used in the numerical model were estimated [8,9]. A summary of these properties is presented in Table 2.  The soil water evaporation curves may be divided into two main stages: 1) a constant evaporation rate during the initial stage, and 2) a decreasing evaporation rate during the second stage. The first stage exhibits a constant slope that is approximately equal to the PE rate. This stage is strongly governed by the atmospheric conditions imposed on the soil-atmosphere interface. During the second stage, the evaporation rate decreases as a function of the soil water retention capacity and soil thickness. The second stage may also be divided into a transition zone and a residual zone. The latter is observed when the moisture content of the soil tends to its residual value, as observed for the samples with thicknesses of 3, 5, and 7 mm. These experimental findings are in general agreement with other studies [6,10].
The mean (μ) wind velocity, air temperature, relative humidity, and irradiance of the seven tests were 1.56 m•s -1 , 41.3 °C, 18%, and 100 W•m -2 , respectively. The standard deviation (σ) values for the wind velocity, air temperature, and relative humidity were 0.09 m•s -1 , 1.3 °C, and 1%, respectively. It is important to note that irradiance was not directly measured but correlated from previous studies using the same climatic chamber and that a constant value of approximately 100 W/m 2 was set for the tests presented herein. For these atmospheric conditions, the mean initial evaporation slope for the soil corresponds to 16.3 g•m -2 •min -1 with a standard deviation of 1.2 g•m -2 •min -1 . The mean value is approximately 0.94 times the PE rate.
For a wide range of atmospheric conditions, a total of twenty-four evaporation tests (AE) were performed on compacted kaolin and twenty evaporation tests were performed on free water (PE). The thickness of the kaolin layers varied between 3 and 20 mm. The initial evaporation rates for all the tests were plotted in a relative humidity (RH) versus PE and wind velocity (V w ) space, as shown in Fig. 4. In this figure, the initial evaporation rates in soil (AE) and water (PE) are assumed to be similar for the similar atmospheric conditions.

Fig. 4. PE rates (results on free water, 20 tests) and initial evaporation rates of compacted kaolin (24 tests).
Equation 1 is the empirical correlation that best fit the initial evaporation rates (i.e., denoted as PE) of the forty-two tests. Where PE is the potential or initial evaporation rate in g•m -2 •min -1 , V w is the wind velocity in m•s -1 , a is a dimensionless fitting parameter calculated as 1.09, e is the Euler's number, and RH is the air relative humidity, in percentage. This expression is valid for RH values greater than 15% and mean wind velocities within the range of 0 to 2.0 m•s -1 , both atmospheric parameters measured near the soil surface (i.e., 5 mm above the soil surface).

Empirical soil evaporation model
The soil evaporation curves obtained in this study are in agreement with the results obtained in other studies [6]. In reference [6], the ratio of Actual Evaporation to Potential Evaporation rates (AE/PE) is presented as a function of the soil moisture availability for Beaver Creek sand tested using a column drying method. Using a similar approach, this study proposes that the AE/PE ratio may be calculated as a function of the mean soil suction for thin soil samples that represent the soilatmosphere boundary. This relationship may be fitted to an inverse sigmoid function, as proposed in Equation 2.
Where e corresponds to the Euler's number, s to the (mean) soil suction, and a, b, and c to fitting parameters based on experimental test results. The parameter a is a function of the suction corresponding to the Air Entry Value (AEV) of the soil. Fig. 5 shows the test results for a 20 mm-thick clay sample (bullets) and the empirical model using the expression in Equation 2 (discontinuous line). The parameters a, b, and c should be obtained from experimental testing. The following fitting parameters were used a = 4000 kPa, b = 0.9, and c = 4.0 in the test used as example (H = 20 mm).

Numerical model for soilatmosphere-structure interaction
A THM numerical model code was written in MATLAB to evaluate soil-atmosphere-structure interactions. The numerical model simultaneously solves a series of partial differential equations (PDEs) of water (i.e., in liquid and vapor phases) and heat flow problems in 3D unsaturated porous media using the explicit method. PDEs for soil stress distribution and stress-strain compatibility at the soil-structure interface are also solved in the numerical model. The structural loading, footing foundation layout, and structural stiffness are considered in the solution.
The comprehensive description of the solution of the numerical model and geotechnical considerations may be found in [8,9]. for different types of geotechnical works. A brief description of the main equilibrium PDEs to be solved is described as follows. It is important to mention that the model proposed in this study evaluates the soil deformation only in the elastic domain and that the effects of soil hysteresis due to drying and wetting cycles are not considered at this stage.

Boundary conditions
The boundary conditions limit the extension of the model and control the transfer mechanisms between the soil-atmosphere and the soil-structure interfaces. The atmospheric conditions and physical laws that govern the exchanges of liquid water, vapor, heat, and loading are imposed on the soil surface. Thus, the soilatmosphere and soil-structure interfaces become the top boundaries of the model. Fluxes along the sides and bottom of the soil model are null. Infiltration rates are imposed from an external database. Evaporation rates are calculated using the empirical expressions proposed in Equation 1 and Equation 2 based on available atmospheric data (i.e., V w and RH) and the SWCC.
The surface of the soil is partially exposed to the atmosphere and partially covered by the building footprint. The building footprint acts as an impermeable barrier, but water vapor may form and condensate beneath the building. Within the limits of the building, the sensible heat is only calculated based on the heat due to thermal emissions and convection and liquid water and vapor flow in the horizontal directions. The soilstructure interface nodes share geotechnical and structural properties. Differences in soil and structure stiffness and differential movements due to soil heave or contraction are calculated, which allows for the calculation of relative displacements between the soil and the foundation.

Heat flux
The sensible heat, q sens , at the soil surface is calculated as the result of the energy balance shown in Equation 3, where q rad is the heat produced by the solar radiation, q th corresponds to the thermal emissions, q conv is the heat due to convection/advection, q rain is the heat due to falling rain, and q evap is the latent heat of evaporation.
The heat flux within the soil mass is calculated by solving Equation 4, where K H is the unsaturated soil thermal conductivity, ρ soil dry density, C H unsaturated soil volumetric heat capacity (e.g., it considers the volumetric heat capacity of the solids and liquid waterair phases), and T is the soil temperature to be solved.

Water and vapor fluxes
The flow of liquid water and vapor through the soil mass is calculated as the volumetric water change for compressible materials, as shown in Equation 5. The solution considers the unsaturated hydraulic conductivity of the soil (k w ), the hydraulic potential (H) of the position (z) of each node, the soil suction (s) based on the Soil Water Characteristic Curve (SWCC), the absolute air (u a ) and vapor pressures (u v ), and the molecular diffusivity of vapor (D v ) through the porous material. The solution of this conservation equation includes changes in the soil degree of saturation (S r ), porosity (n), and temperature (T). Thus, volumetric changes due to shrinkage or swelling may be calculated depending on the infiltration and evaporation balance (i.e., Hydraulic Balance = Infiltration -Evaporation).

Soil-structure interaction
The soil-structure interaction was evaluated by the Winkler's method (1867), in which soil and structure deformations are calculated based on the stiffness of the structural beams (i.e., footings) and the modulus of subgrade reaction of the supporting soil. The basic PDE used to solve the response of the soil and the structure is given in Equation 6. Where EI is the flexural stiffness of the structural element (i.e., footing), K s is the coefficient of subgrade reaction, w is the vertical displacement of the footing, w exp is the vertical displacement of the soil, x is the length of the footing, B is the width of the footing, δ k is a parameter that depends on the soilstructure contact (δ k =1 and δ k = 0, for contact and without contact, respectively), and q is the external load.
The distribution of vertical stresses beneath the foundation system is calculated using the Boussinesq's theory (1885). The displacement (w) at the soil-structure interface is given by the soil deformation due to soil shrinkage or swelling, ε exp , and soil suction, s, curve calculated from the SWCC.
The results of the soil-structure interaction analysis may be combined with the Boscardin and Cording method (1989) [11] to estimate the structural damage. However, this analysis will be presented in further publications.

Results of the numerical simulation
The numerical model was used to evaluate the soilatmosphere-structure interaction for a lightly-loaded, single-story building with a footprint area of 7.5 x 8.0 m 2 . In the model, the soil mass has a volume of 12 x 12 m 2 by 4 m in depth. The building is supported by a continuous footing foundation system as the observed in Fig. 6.   Fig. 6. Results of the numerical model showing the total and relative displacements of the soil and the foundation system.
In this example, the atmospheric conditions correspond to mean values for the city of Cartagena, Colombia (e.g., RH = 75% ±10%, T air = 29 °C ± 5 °C, V wind = 1.5 m/s, Lat: 10° 23' 59''). The model calculates the maximum daily irradiance based on the latitude of the site. However, the irradiance was limited to 50% of the maximum value to account for possible cloudiness. Two dry seasons and two rainy seasons with maximum precipitation values of 1000 mm/year were modelled in 180 days. The loading, dimensions, and bending stiffness (EI) of the structural members were known. The results of the numerical model show the vertical displacements of the building footprint and the surrounding soil (Fig. 6). As expected, sections with higher stiffness exhibit lower relative displacements. Fig. 7 shows the concentration of compressive stresses after a simulation time of 180 days. The clear zones (i.e., white areas) observed near the corners of the building indicate that the contact between the soil and the foundation is only partial. These areas overlap with zones in which the soil deformation due to shrinkage was larger. To compensate for the new vertical distribution of stresses, some sections of the footings take larger loading (e.g., deeper concentration of vertical stresses). Relative displacements between the soil and the foundation and spatial variations of soil saturation degree and soil temperature may also be calculated with the model, but figures were not included.

Discussion
This study simulated the thermal and hydraulic exchanges at the soil-atmosphere interface using atmospheric parameters (i.e., relative humidity, temperature, wind velocity, irradiance) and soil properties (i.e., SWCC, porosity, density) that are easily measured or estimated.
The experimental and numerical approaches proposed in this study for the analysis of the soilatmosphere-structure interaction may be used for different soil types (i.e., fine and granular) and structure geometries (i.e., footings, mats).
Further work is required to validate the experimental and numerical results at larger scales. The elasto-plastic behaviour of the soil may be coupled to the model in future analyses.