Modeling of Evaporation-Driven Multiple Salt Precipitation in Porous Media with a Real Field Application

Soil and groundwater salinization are very important environmental issues of global concern. They threaten mainly the arid and semiarid regions characterized by dry climate conditions and an increase of irrigation practices. Among these regions, the south of Tunisia is considered, on the one hand, to be a salt-affected zone facing a twofold problem: The scarcity of water resources and the degradation of their quality due to the overexploitation of the aquifers for irrigation needs. On the other hand, this Tunisian landform is the only adequate area for planting date palm trees which provide the country with the first and most important exportation product. In order to maintain the existence of these oases and develop the date production, a good understanding of the salinization problem threatening this region, and the ability to predict its distribution and evolution, should not be underestimated. The work presented in this paper deals with the Oasis of Segdoud in southern Tunisia, with the objective of modeling the evaporation-driven salt precipitation processes at the soil profile scale and under real climatic conditions. The model used is based on the one developed and presented in a previous work. In order to fulfil the real field conditions, a further extension of the geochemical system of the existing model was required. The precipitated salts considered in this work were halite (NaCl), gypsum (CaSO4) and thenardite (Na2SO4). The extended model reproduces very well the same tendencies of the physico-chemical processes of the natural system in terms of the spatio-temporal distribution and evolution of the evaporation and multiple-salt precipitation. It sheds new lights on the simulation of sequences of salt precipitation in arid regions. The simulation results provide an analysis of the influence of salt precipitation on hydrodynamic properties of the porous medium (porosity and permeability). Moreover, the sensitivity analysis done here reveals the influence of the water table level on the evaporation rate.


Introduction
Saline water evaporation from porous media is involved in several applications such as building industries, agriculture and water resources management [1][2][3][4][5].
The evaporation of saline water from soils results in the precipitation of the dissolved salts, leading, therefore, to groundwater and soil salinization. This environmental problem is more pronounced in arid and semiarid regions where, apart from dry climatic conditions, the irrigated water, often salt-riched and of poor quality, provides an additional source of salts to the soil. with the developed model for a soil profile typical of the oasis parcel and using realistic climatic data provided by the meteorological station corresponding to the region subject of this work. In what follows, the formulation of the mathematical and geochemical models is detailed. The simulation results related to the model application are discussed in the last part of this paper.

Model Assumptions and Conceptualization
The processes of evaporation driven salt precipitation in porous media are influenced by the initial chemical composition of the evaporating solution leading to favorize the precipitation or dissolution of one salt type in disregard of the other. Hence, to be able to describe the real field system, a good understanding of the brine chemistry is required for the model conceptualization.
The chemical composition of the brine solution is given in Table 1. It results from the equilibrium between irrigated and drainage waters chemistry. This equilibrium indicates that HCO 3− ions are eliminated in drainage water. Moreover, reference [40] indicated that the irrigation water used for Segdoud oasis is very deep water which brine is very poor in bicarbonate ions. Hence, our first assumption in the model concept is to avoid HCO 3− in the chemical system. As can be seen from Table 2, the Segdoud oasis groundwater is Na + − Ca 2+ − Cl − − SO 2− 4 enriched. The concentrations of K + and Mg 2+ are very low compared to the other ions. For that reason, our second assumption is to ignore the presence of K + and Mg 2+ and to consider a four-ions based chemistry with a correction of the electro-neutrality of the solution. Table 1. Chemical composition (in mg/L) of Segdoud brine (March 1995) taken from [42]. Based on this assumption, we did preliminar batch reaction calculations with the geochemical model PHREEQC [43]. These calculations are helpful for problems where the possibility of mineral dissolution or precipitation needs to be known and analyzed, which is the case in this work. From a brine solution analysis, PHREEQC calculates the distribution of the ionic species and provides saturations indices for minerals indicating whether a mineral should precipitate or dissolve. Table 2. List of minerals prior to precipitate after evaporation of Segdoud brine.

Minerals Chemical Formula
Calcium Sulphate (Gypsum) In the case of Segdoud brine which composition is given in Table 1, different minerals can be formed. These are listed in Table 2. The sequence of precipitation of these salts by evaporating the initial solution is shown in Figure 1. It depicts the evolution of the molality of the ionic species as well as the formed minerals as a function of the concentration factor of the initial solution during evaporation.
According to the batch reaction calculation results shown in Figure 1, it is shown that only three minerals (salts) should precipitate during the evaporation of the initial brine characterizing the oasis of Segdoud. These salts are gypsum, thenardite and halite, respectively. Therefore, the model concept developed for the modeling of real evaporation driven salt precipitation from the soil profile representing the oasis of Segdoud should take into account the simultaneous precipitation of these salts. The model concept developed in this aim is described by Figure 2.
. There is exchange of dissolved ions between liquid and solid phase through precipitation and formation of Gypsum (G), Halite (H) and Thenardite (T).

Geochemical Modeling
The chemical system considered for the developed model concept accounts for the precipitation of gypsum, halite and thenardite. The developed model concept is an extension of the model developed by E.Mejri and described in [41] with the required components and equations. Hence, we give here the description of the geochemical processes related to the precipitation of thenardite. All details describing the kinetics of gypsum and halite formation can be found in [41].
Thenardite (anhydrous) is one of two widely known forms of sodium sulphate, the other one being mirabilite (Na 2 SO 4 · 2H 2 O). According to the literature [19,[44][45][46], sodium sulphate is considered one of the most harmful salts causing the degradation of porous building materials. It is widely used in laboratory experiments dealing with rock weathering problems. Thenardite forms according to the following chemical reaction: with an equilibrium constant K eq,T given as: where [i] is the activity of the ionic species i, calculated using the Truesdell and Jones model [41]. At 25 • C and constant ph, the solubility of thenardite in water is 160 g/L. The chemical processes are described following the kinetic approach. This suggests that thenardite forms with a precipitation rate r precip,T calculated as follows: The parameters of this equation are found in [41]. The key parameter in this equation which indicates if thenardite should precipitate or dissolve in the liquid solution is the saturation index Ω T . It is calculated as: where the solubility product of thenardite is equal to K sp,T = 3.764 [ mol 2 Kg 2 ]. Combined with the evaporation dynamics, the analysis of this multiple salt precipitation system is of great interest because the considered salts have different characteristics (anhydrous and not, solubility, ...) and the dissolved ions participate in more than one chemical reaction. This would create a competition between the formation of the different salts. These processes are discussed later on with the help of the numerical simulations.

Mathematical Modeling
The reactive transport model developed in this paper is an extension of the model developed in [41]. While the same existing transport equations of the components were reused here, an additional equation was needed to describe the transport of thenardite. It is given as follows: where φ T corresponds to the volume occupied by the precipitated thenardite (solidity). Moreover, the source/sink terms of the transport equations describing the coupling between transport and chemical processes must be modified taking into consideration the common-ion salt precipitation. Table 3 summarizes the source/sink terms of the transport equations for each component of the multiphase-multicomponent system, where r precip,G and r precip,H are defined in [41]. The ionic species dissolved in the brine solution are subject to transport in the porous media and during the course of evaporation, their concentration increases leading after some time to the formation of salt crystals within the void space of the porous medium. The decrease of the initial porosity of the modeled porous structure is evaluated with the following equation: where φ m is the solidity of the precipitated mineral m, m ∈ {G, T and H}. Besides, the same equation accounting for the reduction of the permeability of the porous medium resulting from the reduction of its porosity is used. The solidity of thenardite φ T is added to the previously considered primary variables to solve the non-isothermal multicomponent reactive transport model. No modifications were done on the numerical methods used to discretize the system of equations [41]. Table 3. Source/sink terms of the transport equations of each component.

Component Source/Sink Term
Water

Location
The study area is the Oasis of Segdoud located in the southwest of Tunisia (34 • 14 N; 8 • 10 E), 30 Km north-east of Tozeur and 70 Km south-west of Gafsa (see Figure 3). Segdoud is limited to the north by the Negueb and Chouabine mountains and to the south by chott el Gharsa. The Oasis of Segdoud is a modern oasis which has been created in 1989 by the country policies in order to enhance the economical development of the country by maximizing the production and exportation of dates. The oasis covers a surface of 160 ha and is subdivided into 107 plots, each covers an area of 1.5 ha. The oases of Segdoud are continental with a mean planting space of 10 × 10 m [47]. The date palm is the principal vegetation in these oases with the preference of Deglet el Nour variety.

Climatic Conditions
The oasis of Segdoud belongs to the bioclimatic saharan zone. The weather is characterized by an annual average temperature of 20 • C. The warmest month is July with a mean temperature of 30.6 • C while the coolest month is January, with a mean temperature of 10.6 • C. We plot in Figure 4 the evolution of the mean temperature during the year 2012. The data are obtained from the National Institute of Meteorology and collected at a daily time step from the meteorological station located in Tozeur.
In Tozeur, the precipitations are low and irregular with an annual amount of 96 mm on a total of 29 days. March is the month with the most precipitations with a mean value of 20.3 mm per year. The month with the least precipitation is July with an average of 0 mm yearly. The weather in Tozeur is hot and dry during summer and is characterized by a high evaporation rate. The annual reference evapotranspiration in Tozeur is equal to about 1643 mm. Figure 5 depicts the evolution of the evapotranspiration in Tozeur during 2012. These calculations are done with the refererence evapotranspiration calculator ET0 using the FAO (Food and Agriculture Organization) Penman-Monteith Equation [48]. They are based only on weather data and depend neither on crop characteristics nor on soil parameters. The FAO Penman-Monteith equation is given as follows:

Chemical Composition
The soils in the oasis of Segdoud are 96% of sandy type. The variability of the soil texture is very low from one soil layer to the other [42]. The saturated hydraulic conductivity of the oasis soil varies between 0.2 and 30 cm · h −1 and its measured electric conductivity lies between 8 and 30 dS · m −1 . Such a value allows to classify the oasis soil within the saline soils [42].
Hydrographically, the oasis of Segdoud is connected to Chott el Gharsa, which contains three superposed aquifers: The Complex Terminal (CT), the Continental Intercalaire (CI) and the shallow groundwater aquifer.
Based on experimental investigations, reference [42] showed that the depth of the shallow groundwater table varies from 0.9 to 1.8 m. Water flows in the oasis from north to south and from east to west with a mean infiltration rate of 2.2 m/day and a hydraulic gradient lying between 0.7% and 1.5%. The salinity of the oasis groundwater is high and results from an equilibrium between the salinity supply from irrigation water and salinity removing by the drainage and groundwater lateral flow. These waters are rich on different ionic species, namely Na + , Ca 2+ , Mg 2+ , K + , Cl − , SO 2− 4 and HCO 3− [42].

Management of the Oasis of Segdoud
The water management in the oasis of Segdoud started in 1989 after its creation. At that time, the effective irrigated area was about 10 ha from the total area of 160 ha. This irrigated surface has increased to reach 20 ha in May 1995 [38].
The irrigation method is a traditional system of pipes, concretes and earthen canals. The irrigation water is of poor quality and often very saline, as reported by Askri [36]. In Table 4, the chemical composition of water used for irrigation in October 1995 is given. Irrigation water is pumped from two deep wells Segdoud CT2 and Segdoud CT3 and is conducted to the basins in which the palm trees are [36]. The total water inflow by irrigation is 76 L/s. The duration of irrigation is 3.3 h of water per hectare with a frequency which depends on the season: Every 7 days in spring and summer and every 15 days in autumn and winter [36].
Near to the irrigation system, a drainage system was also installed in the oasis in order to discharge the surplus of saline water from the water table [38].

Problem Setup
In the work of [36], it has been reported that, in the oasis of Segdoud different processes take place at the oasis parcel scale. These processes consist in evapotranspiration, rainfall, irrigation, infiltration, percolation to groundwater, groundwater lateral flow and drainage. In our numerical experiment, we do not consider, however, all these processes to avoid complex calculations and long calculation time. The first focus point of our analysis is the evaporative multiple salt precipitation. The other processes will be considered in future investigations.
To do so, we considered a two-dimentional homogeneous porous column of 1.1 m of height and 0.5 m of width and meshed with a 41 × 20 grid, as schown in Figure 6, which also depicts the simulation setup and the initial and boundary conditions. The properties of the porous medium are typical of the Segdoud oasis soil. These are given in Table 5. The bottom and sides of the column were closed. Within the porous medium, a water table level is set at 0.5 m, above which the saturation in the domain corresponds to a hydrostatic pressure distribution. The chemical composition of the initial brine solution is the one given in Table 2, but converted into mass fractions.  At the top of the porous medium, an artificial domain was introduced in which we defined the time dependent boundary conditions. These are given for the temperature in order to calculate the evaporation from the soil profile. The values of temperature are hourly data related to the month of July from the year 2012 with the most complete climatic data. This choice is explained by the fact that July is the warmest month during which the evaporation is the most important in the year. These values are plotted in Figure 7 and introduced as boundary condition at hourly calculation time steps. The first value corresponds to 00:00 AM and the fluctuations are explained by diurnal variations of the temperature.

Simulation Results and Discussion
We discuss in this section the results of the numerical experiment carried out with the previously described setup. The results shown here are limited to 10 days simulation time. In Figure 8, the spatial distribution of the water saturation in the soil profile for different simulation times is given. The dashed line separates the porous medium from the artificial domain.
As it can be clearly seen in this figure, the water level decreases progressively and homogeneously in the unsaturated zone, which is explained by the evaporation from the porous medium driven by the temperature and vapor gradients. A longer simulation time would be needed for the drying front to reach the water table level set in the middle of the porous column. As long as the porous column dries out, the liquid solution is transported from deeper level to the evaporating surface. Influenced by the decrease of water content, the concentration of ions in the liquid moisture increases consequently. This can be clearly seen in Figure 9 which illustrates the temporal evolution of the ions concentration along a vertical cross section of the domain. Analyzing this figure, we see that the concentration of all ions is the highest at the top of the porous column, where water evaporates, compared to the lower part of the domain. Moreover, it is observed that the concentration of Na + and Ca 2+ decreases between 3 and 10 simulation days. This decrease is explained by the supersaturation of the solution with respect to the solubility concentration of the dissolved minerals indicating the beginning of their formation. Since Na + is a common ion in two precipitation reactions, the decrease of its concentration will be the most notified. Furthermore, the initial concentrations of Cl − and SO 2− 4 are higher than for Na + and Ca 2+ which explains that the solution remains concentrated in these two ions despite the beginning of salts precipitation. Figures 10-12 describe the spatial distribution of the volume fraction of the formed gypsum, thenardite and halite, respectively, for different simulation times. Here, it is important to mention that the initial brine solution was oversaturated in dissolved gypsum with an initial saturation index Ω G = 1.12. This results in the presence of precipitated gypsum already everywhere in the domain even before drying starts. After some time, the volume fraction of the precipitated gypsum at the top of the domain becomes lower compared to the rest of the domain. This is because of the simultaneous formation of thenardite and therefore the consumption of SO 2− 4 for that reaction. As for thenardite and halite, they start forming at the evaporating surface where the ions concentration is the highest. This is driven by the drying of the porous column. Since the evaporation driven multiple salt precipitation mainly takes place at the top of the domain, we focus in what follows on the analysis of the processes at a chosen point situated in the middle of the evaporating surface, since the processes are homogeneous horizontally.  It is important to mention, that in Figures 10 and 11, no salt fingering is observed in comparison to what has been reported in other presious studies [49][50][51]. This is explained by the absence of density effect considerations in the developed numerical model. This will, however, be taken into account in future work, where density variation related to salt concentrations change will be involved in the model.
In Figure 13, we plotted the time evolution of the precipitated volume fraction of the given salts at that representation point. We clearly see, that gypsum started forming from t = 0, as explained previously. The gypsum formation is followed after a while by halite and thenardite crystallization.  The results provided by this analysis are not, however, comparable to those calculated using PHREEQC calculator in terms of the time of apparition of the first crystals of halite and thenardite. In fact, the batch reaction calculations performed with PHREEQC are based on the equilibrium approach and not the kinetic one used in our model. This latter supposes that, the higher the precipitation rate constant of a mineral is, the faster would this mineral precipitate under drying conditions. In this numerical example, we used kinetic rate constants for halite and gypsum of the same order of magnitude. Hence, the simultaneous beginning of their crystallization was expected.
Salts precipitation has a direct consequence in clogging the void space of the porous material and, therefore, the reduction of its initial porosity. Such a result is clearly seen in Figure 14. This figure shows the temporal evolution of the initial porosity at the same representation point used previously. The sharp decrease of the porosity occuring at t * is explained not only by the beginning of formation of halite and thenardite, but also due to the higher volume fraction of halite (reaching 2 × 10 −5 ) compared to gypsum (3 × 10 −6 ). In Figure 15, the temporal evolution of the real evaporation rate simulated with the developed model is given and compared to the evaporation rate profile obtained from the work of [29]. Initially, the evaporation rate increases due to the temperature and water vapor gradients between the porous medium and the artificial domain. This is similar to the first stage SS1 of evaporation of saline water. The transition to stage SS2 is marked by the increase of the evaporation rate. Two curve slopes are observed: A first sharp decrease, called here D1 followed by a second lower decrease, called D2. At t = 0.5 day, the evaporation rate drops. From the one hand, at that time, the porosity decrease is not high enough to cause the decrease of the evaporation rate. From the other hand, as long as the porous column dries out, the ions concentration increases. This results in the increase of the osmotic potential of the solution, which in turn leads to the decrease of the saturated vapor pressure and therefore the potential of water molecules to move from liquid to gas phase. Hence, the increase of ions concentration in the brine solution is high enough to cause the sharp decrease of the evaporation rate during D1.
The second decrease D2 is less important and caused by the clogging of the pores filled with the precipitated salts, which causes the reduction of the evaporating surface. From the end time of the second decrease D2 on, the evaporation rate is maintained low and constant. This indicates that the water vapor flow is diffusion driven. This result matches pretty well with the profile given in the upper corner of Figure 15 [29].
Furthermore, the evaporation rate obtained from these simulations is compared with the calculations done with the calculator ET0 using the same initial climatic conditions (in terms of temperature and humidity). This comparison is given in Figure 16. We clearly notice from this figure, that at the beginning the results are in good agreement with each other. At that time, the water saturation at the top of the simulated porous medium is high and the evaporation is driven only by the temperature and vapor pressure gradient. From the second day on, water saturation decreases and salt concentration increases leading to the decrease of the simulated evaporation rate while the evaporation rate calculated by ET0 is maintained high. So far, the developed model enables a good description of the evaporation driven mutiple salt precipitation processes taking place in the soil profile characterizing the oasis of Segdoud. Apart from that, the depth of the water table and its seasonnal fluctuations play an important role in the soil salinization of the oasis of Segdoud, as reported in [36]. Hence, we performed a sensitivity analysis on this parameter to provide a better understanding of the influence of the water table level on the evaporative salt precipitation inducing soil salinization. We performed in this aim numerical examples with two water table levels: L 2 = 0.4 m and L 3 = 0.3 m as shown in Figure 17. In Figure 18, we plot the simulation results in terms of the temporal evolution of the evaporation rate for the given three water table levels (level L 1 = 0.5 m corresponds to the first simulation test). Based on these plots, we observe that the evaporation rate is the highest for the higher water table level. Actually, the nearest is the water table to the evaporating surface, the more water is available and exposed to atmospheric conditions. The evaporation rate is maintained the longuest time constant for the deepest water table level. This result is supported by the results depicted in Figure 19 which describes the evolution of the initial porosity of the porous medium at the evaporating surface for the tested water table levels. We see that the porosity decrease becomes important as the water table gets closer to the evaporating surface. This indicates that soil precipitation is more important when the water table is not deep, which confirms the results given by [42].

Summary and Outlook
In this paper, we presented the model developed to simulate and analyze evaporation driven salt precipitation processes at the soil profile scale and under real climatic conditions related to the oasis field of Segdoud (Southern Tunisia). The developed model is based on the existing model [41] and has been extended to simulate the dynamics and kinetics of drying and mixed salt precipitation processes.
The issue of evaporation salt precipitation cases where more than one salt is precipitating, is lack of treatement in the literature [7,26,29]. This work is, hence, considered a novelty in this direction since it provides a good description of processes dealing with evaporative multiple salt precipitation in porous media. In fact, Several simulations have been performed with the developed model. In all simulations, far field boundary conditions in temperature were used and taken from the National Institute of Meteorology (INM, Tunis). The model reproduces very good the different processes related to the soil evaporative salinization of the oasis of Segdoud.
Besides, the simulation results given in terms of the evaporation rate curve and the effects of water table leven on the evaporated water, are in good agreement with results published in [33,37,42,52].
So far, the developed model concept is restricted to the porous media at the top of which an artificial domain was introduced to define boundary conditions of temperature and water vapor. Many works existing in the literature, demonstrated, however, that the evaporation processes are highly influenced by different weather conditions. These are e.g wind velocity, turbulence, relative humidity and radiation. Therefore, a coupling of the developed porous medium model with a free-flow model has to be considered to describe the evaporation processes with more accuracy.
Furthermore, for the sake of simplicity, the implemented geochemical model accounts for the simultaneous precipitation of three salts. In reality, there exist cycles of evaporation leading to dissolution and recrystallization of the salts resulting in more complex precipitation phenomena and competition between the precipitated salts. The precipitated salts themselves, can react with the solid matrix causing its dissolution and fracturing. These phenomena have not been taken into account in this work but need to be investigated in future research. Moreover, the processes investigated in this paper are density dependant because of the variations in water saturation and dissolved minerals concentration. Hence, it would be interesting to point out in future work the processes of density recirculation in the simulated soil profile.