Simulating the drying behavior of rammed earth columns

Rammed earth construction provides an efficient alternative construction material to limit energy consumption and CO2 emission. It possesses various characteristics of a sustainable material, but its mechanical capacity is sensitive to humidity variations. It is, therefore, important to better understand water transport within rammed earth when subjected to varying ambient conditions. In this context, the present work aims to analyze the hydraulic behavior of a reduced piece of rammed earth wall consisting of a column of size 14cm x 14cm x 30cm subjected to drying from the initial compaction water content in an indoor environment. The columns had transient non-uniform relative humidity, which was measured in-situ. Thermo-hydraulic coupled numerical modeling was developed using realistic atmospheric boundary conditions. The material and hydric parameters were chosen from an experimental study previously performed at the material scale. A parametric study was performed in order to evaluate the sensitivity of the modeling to both material parameters and boundary conditions. The results of the numerical simulations were highly sensitive to parameter values used for the water retention curve and the surface mass transfer coefficient, with a satisfactory matching of experimental results only achieved after adjusting initial estimates of relevant parameter values.


Introduction
Rammed earth is a sustainable construction material that can help to reduce the carbon footprint and limit the energy consumption in the construction engineering domain. It has lower embodied energy than concrete or steel [1], it is recyclable in its unstabilized state and has desirable hygro-thermal properties [2] [3]. In spite of these advantages, it is vulnerable to water ingress which affects the strength and durability. After its preparation at optimum moisture content, rammed Earth dries, leading to a rise in strength. This gain in strength is attributed to additional suction-induced cohesion in the unsaturated state [4][5] [6]. Conversely, when it is subjected to wetting (e.g., rainfall or inundation), a reduction in strength is observed. During the lifetime, rammed earth structure is subjected to changing humidity conditions (affecting the matric suction), which induces transient moisture transfers through the wall. Thus, in order to predict the mechanical response, it is important to estimate the hydraulic conditions of the material accurately. François et al. 2017 [7], performed hydraulic modeling of rammed earth walls by imposing a water pressure at the wall faces (and thus suction) as the boundary condition. With an applied suction as a boundary condition, the surface of the wall reaches the final suction state immediately and is supposed to be in equilibrium with the surrounding atmosphere. However, equilibrium condition between external conditions and the material surface is reached over an infinite time asymptotically. * Corresponding author, e-mail: parul.chauhan@univ-smb.fr Thus, in the present work, a coupled thermo-hydraulic finite element modeling is done using a more realistic atmospheric boundary condition. This contribution aims to simulate the progressive drying internal of a piece of rammed earth wall subjected to drying condition, from its initial compaction water content. The atmospheric boundary condition is expressed as a flux boundary condition. Firstly, details about the material properties are presented. Experimental methods and results for 1-D drying of the column are explained. A parametric study is performed to analyze the sensitivity of the modeling to both material parameters and boundary conditions. A comparison between the experimental and numerical results for the four parameter sets is described.

Material
The material studied here comes from a rammed earth construction site in Auvergne-Rhone Alpes region of France. The particle size distribution curve was plotted using wet sieving and hydrometer analysis according to French standard NF P 94-057 [8]. The PSD curve in figure 1 shows that it contains 7% gravel, 30% sand, 49% silt, and 14% clay. Its index properties are: liquid limit w l = 27.4%, plastic limit w p = 16.4%, and plasticity index I p = w l − w p = 11.0%. It is classified as low plastic silt (I p < 12%) according to the French classification GTR (Guide de Terrassements Routier) for fine soils. The activity of clay (A c = I p / f ) defined as the ratio of plasticity index (I p ) and percentage of soil passing 2 µm sieve (f) was equal to 0.78. The high value of Specific surface area (S sp = 14.7 m 2 /g) and Cation exchange capacity (CEC = 2.6 cmol/g) suggests a very low percentage of swelling clays. Figure 1. Particle size distribution of the soil used for sample manufacturing [9] Optimum moisture content (OMC) of 12.5 % was obtained from the Standard Proctor test. The soil was mixed at this water content and left for water content equalization in a sealed container for 24h at least.

Experimental study on columns subjected to drying
One column of size 14cm x 14cm x 30cm was manufactured by dynamical compaction of soil prepared at OMC in between removable wooden frameworks (figure 2b). The soil was poured loosely in four layers of about 10-15cm into the framework and rammed by moving the hammer from the edge to the center of the specimen until reaching a homogeneous and leveled surface. After compaction, the thickness of each layer was typically around 6 to 10cm. When the required height was achieved, the upper part of the mold was removed, and the surface was leveled. The sample had an average dry density (ρ d ) of 1907 kg/m 3 . The range of dry density for earthen structures is 1700 kg/m 3 to 2200 kg/m 3 [10] and the obtained dry density lies within this range. The prepared samples were covered with paraffin from the top, bottom, and two parallelfaces imposing a unidirectional drying. The paraffin was put in that manner to reproduce the boundary conditions of a representative elementary volume taken from a bigger rammed earth wall (figure 2a). Three capacitive SHT sensors for the measurement of relative humidity and temperature were added in the middle of the third layer. The SHT sensors were placed at a distance of 3cm, 5cm, and 7cm from a drying face (figure 2c). The sample was left to dry in an indoor environment. The values of ambient relative humidity and temperature were measured at regular intervals until day 61. The mean value of ambient RH and temperature recorded were 65 % and 15.5 • C respectively. The transient profiles of relative humidity within the sample are plotted in figure 3. It appears that the relative humidity at 3cm from the surface was not affected until day 13. Similarly, at a distance of 5cm and 7cm (center of the column), no change in RH is observed until day 24 and 31 respectively. After reaching this stage, a significant decrease at a substantial rate is observed in the RH, until reaching a final stage of RH = 73 %, 74 %, and 78.5 % at 3cm, 5cm, and 7cm respectively at day 61.

Numerical study
Numerical simulations of the drying process of rammed earth were performed using CODE_BRIGHT (Olivella et al. 1996 [11]), which is finite element code to carry out thermo-hydro-mechanical (THM) coupled problems in geological media. Here, Thermo-hydraulic analysis was performed considering a rigid solid phase and an immobile gas phase. Thus, displacements (u) were considered to be 0 and gas pressure (P g ) was assumed to be atmospheric (P g = 0.1 MPa). Energy balance (solving for unknown temperature, T) was taken into account. Vapor diffusion into the gas phase was considered, and in fact, diffusion has a significant effect on drying. The geometrical model (figure 2d) considered was a 2D column of height 0.3m and a width of 0.14m to represent the column studied. The mesh was refined only along the direction of drying (x) since it is a 1D problem. No mesh elements were required in the y-direction. However, CODE_BRIGHT requires a minimum of 3 elements in each direction. A mesh of 50 x 3 (150) quadrilateral (Q4) elements was used. The mesh generated had a '0.1' concentration of mesh close to the edges in the x-direction. The total number of nodes was 204.

Theoretical aspects of CODE_BRIGHT
The retention behavior of rammed earth was defined by using the Van Genuchten 1980 [12] (VG) model. The Van Genuchten expression for effective degree of liquid saturation (S el ) is defined as: where, S l is the actual degree of liquid saturation, α(MPa) −1 , n, and m are model parameters, s is matric suction(MPa). Here, the maximum (S ml ) and residual (S rl ) degree of saturation were considered as 1 and 0 respectively, so S el = S l . m is related to n as m = (n − 1)/n.
The unsaturated hydraulic conductivity function was defined as k = k r k sat , where k r is the relative hydraulic conductivity and k sat is the saturated hydraulic conductivity. Mualem hydraulic conductivity model (1976) [13] was used to express k r as a function of S r coupled with Van Genuchten retention model, which gives: The intrinsic permeability (k i ) can be related to the saturated hydraulic conductivity (k sat ) by the expression: where, ρ l is the density of water and µ w is the dynamic viscosity of water.  The convection equation for the liquid flow rate (q l in m/s) used was the Darcy's law: where, µ l (MPa.s) is the liquid viscosity, P l is liquid pressure (MPa). The diffusive flux of water vapor in gas phase (i w g ) is evaluated using Fick's law of vapor diffusion: where, τ= tortuosity, φ =porosity, ρ g = gas density, S g = 1 − S l is the gas degree of saturation, D w g (m 2 /s) is diffusion coefficient of water vapor in gas phase, ω w g (kg/kg) is the mass fraction of water vapor in gas phase, D(m 2 /s) and d are the parameters of the model. The conductive heat flux i c (W/m 2 ) is computed using Fourier's law: is the thermal conductivity, λ solid , λ gas , λ liq Wm −1 K −1 are the thermal conductivity of solid, gas and liquid phases respectively.

Material parameters
The data points for the soil water retention curve were taken from a previous study on this soil [14]. The VG model parameters were evaluated as α = 14.1 MPa −1 , n=1.287, and m=0.223 and the retention curve is shown in figure 4 labeled as VG. This retention curve is not accurate since there is no data for s<3.8MPa, therefore, an additional fictitious data point was considered at the 'ascompacted state' for which we determined S l = 0.805, and we suppose a value of suction of 0.5 MPa (much lower than 3.8 MPa, value for which S l is equal to 0.30). This was done to check the sensitivity of the modeling to the retention behavior. The parameters considering the additional data point were evaluated as α = 2.763 MPa −1 , n=1.4, and m=0.2857. It is labeled as VG modified (figure 4). The saturated hydraulic conductivity was equal to 3.275 * 10 −9 m/s, and the intrinsic permeability was evaluated using equation 3 as 2.9147 * 10 −16 m 2 . Due to uncertainties in different material parameters, it was decided to conduct a parametric study, which will be explained in the further section. The parameters used for the diffusive flux of water vapor and the conductive heat flux are presented in table 1.

Initial conditions and boundary conditions
From the experimental data, the initial degree of saturation after compaction at OMC was 0.805. The initial condition for liquid pressure (P li ) of the numerical model was chosen according to the retention curve adopted. For the retention model VG, the initial suction corresponding to S l = 0.805 was s= 0.105 MPa, as the gas pressure is P g =0.1MPa, the liquid pressure (P li ) applied to the model was equal to -0.005 MPa. Similarly, for the VG modified retention curve, at S l =0.805 the suction is s = 0.396 MPa which corresponds to P li = −0.296 MPa. The initial temperature was T=15.5 • C. When an imposed liquid pressure (and hence suction) is applied as the boundary condition, the boundary immediately reaches the final suction state and is supposed to be in equilibrium with the atmosphere. This condition is not realistic since the equilibrium at the surface is approached asymptotically. Thus, in this study, a more realistic atmospheric boundary condition is used. It allows imposing boundary conditions in terms of evaporation and heat exchanges. Flux boundary condition was used to express these phenomena for air and energy in terms of the state variable P l and T respectively, and the ambient conditions such as relative humidity, atmospheric gas pressure, temperature, and air velocity. Here the data related to solar radiation and rainfall have not been used to use this boundary condition for an indoor environment. In order to simulate natural convection, the flux due to evaporation has been modified since this atmospheric boundary condition is generally used for an outdoor environment. In CODE_BRIGHT, the evaporation flux E (kgs −1 m −2 ) is defined by an aerodynamic diffusion relation: The evaporation flux can also be defined as [15] where, ρ va and ρ v are the absolute humidity of the atmosphere and at the boundary respectively, k is the von Karman's constant, s f is a stability factor, v a the wind velocity,

Parametric study
As mentioned in the previous sections, there is uncertainty in the values of both material parameters and boundary conditions. Thus, a parametric study was done to understand the influence of a single parameter or combination of parameters. Firstly, a parameter set (PS1) was chosen, considering the default values of all the parameters. VG retention curve and value of h c =0.85Wm −2 K −1 as suggested by Kusuda 1983 [16] were selected. In parameter set PS2, VG modified was selected to understand the effect of change in retention behavior on hydraulic modeling. Other parameters were kept constant. In parameter set PS3, a higher value of h c =25Wm −2 K −1 was used to study the effect of increased h m keeping the other parameters constant. Finally, in the parameter set PS4, VG modified, and h c =25Wm −2 K −1 was used in combination. It is to be noted that the parameters for vapor diffusion and conductive heat flux are the same in all the simulations (table1), and a value of k i used was equal to 2.9147 * 10 −16 m 2 . These four parameter sets are mentioned in

Results and discussion
In this section, the results of the four simulations based on the four parameter sets are presented. Figure 5 and 6 show the isochrones of suction and degree of saturation respectively across the width of the rammed earth column for PS1. From the suction profiles, it can be seen that suction at the boundary does not reach to the final equilibrium state immediately and is gradually increasing over time. Thus, there are lower gradients of suction between the boundary and inside of sample than would occur if the final suction was applied to the boundary from t=0. From figure 5, it appears that even after 31 days of the drying period, the suction across the width is hardly affected, but from the degree of saturation profiles, it is clear that a significant amount of water has been taken out. S l at the boundary decreases from 0.805 to 0.49, whereas at the center of the sample, it decreased to 0.65. Also, after 31 days, suction at the boundary and the middle is equal to 858.2 kPa and 279.8 kPa respectively compared to initial and final values of 105 kPa and 55.85 MPa, respectively.  Figure 7 shows the comparison of the experimental and the numerical results for the parameter set PS1. Up to 70 days, the numerical simulations predict little change of RH at the three sensor locations, but as discussed before, a significant reduction of S l has already been predicted. This predicted behavior of RH was observed because of the non-linear relation of both S l − s and s − RH. It can be seen that the simulated drying process is too slow compared to the drying measured: 200 days are needed instead of 61 days. Although the drying process takes a longer time, the pattern of drying is similar. The difference of RH at 5cm and 7cm is lower compared to 3cm and 5cm. This gradient across the width first increases and then decreases as the drying process approaches its final stage. Figure 8 shows the predicted RH variation for parameter set PS2, where the retention behavior is changed while the atmospheric boundary condition applied is the same. The drying process was relatively fast, and the gradients of RH across the width were lower compared to PS1. Considerable changes in RH were observed at the beginning of the simulations. From 70 days onward, steeper rates of RH were observed across the width at 3cm, 5cm, and 7cm. Still, there is a significant difference in the simulation results compared to experimental results. It is to be noted that due to change in the retention behavior, the initial conditions of simulations were also changed. Drying was started from a higher value of initial suction. This had a significant impact on the drying process.
In the parameter set PS3, the retention behavior is unchanged while the atmospheric boundary conditions are changed, and a higher value of surface mass transfer coefficient is applied. Figure 9 shows that the pattern of drying is similar to PS1 but is accelerated, which seems to be closer to the reality in comparison to PS1 simulation. The whole width of the column reached close to the final value of RH in 200 days. There were significant gradients of RH across the width, which is similar to PS1. This shows that the higher gradients were due to the retention curve as for PS2, lower gradients were observed. In order to understand the combined effect of change in retention properties and higher surface mass transfer coefficient, both the material property and boundary condition was changed in the parameter set PS4. Figure 10 shows that significant changes in RH occurred from the beginning of the simulation across the width. The results at the end of 61 days were closer to the experimental results, although the initial behavior of drying was not reproduced accurately. The RH for the whole width reached the final value of 65% in less than 150 days. The suction values at the boundary in the beginning stages of the simulation ascended faster compared to the three simulations before. The suction at the boundary was as follows: at t=1 day, s=16.14 MPa, at t=2 days, s=36.3 MPa and at t=5 days, s =47.81 MPa. For this reason, very high gradients of suction were observed across the width. Similarly, the degree of saturation at the boundary was close to its final value in the first few days of the simulation. It is clear that the simulations are very sensitive to both the retention curve and boundary conditions.

Conclusions and future work
In the present study, the experimental method and results of an element of the wall subjected to drying were presented. After a certain period of time, a significant decrease in relative humidity at a substantial rate was observed. A coupled thermo-hydraulic modeling using CODE_BRIGHT was done. Atmospheric boundary conditions were used due to which the suction at the boundary reached its final stage gradually, which is more realistic than an imposed suction condition at the boundary. The parametric study performed shows that the results of the hydraulic analysis were very sensitive to the soil water retention behavior and the surface mass transfer coefficient, i.e., boundary condition. Using the modified retention curve, higher intrinsic permeability, and higher mass transfer coefficient, it was possible to reproduce the experimental results after 61 days of drying fairly, but the fitting to initial phase of drying period was inadequate. The uncertainty in the retention behavior can be eliminated by having additional data points for suction less than 3.8 MPa. Also, the value of the surface mass transfer coefficient can be calibrated experimentally.