Modelling Shallow Groundwater Evaporation Rates from a Large Tank Experiment

A large tank (1.4 m x 4.0 m x 1.3 m) filled with medium-coarse sand was employed to measure evaporation rates from shallow groundwater at controlled laboratory conditions, to determine drivers and mechanisms. To monitor the groundwater level drawdown 12 piezometers were installed in a semi regular grid and equipped with high precision water level, temperature, and electrical conductivity (EC) probes. In each piezometer, 6 micro sampling ports were installed every 10 cm to capture vertical salinity gradients. Moreover, the soil water content, temperature and EC were measured in the unsaturated zone using TDR probes placed at 5, 20 and 40 cm depth. The monitoring started in February 2020 and lasted for 4 months until the groundwater drawdown became residual. To model the groundwater heads, temperature, and salinity variations SEAWAT 4.0 was employed. The calibrated model was then used to obtain the unknown parameters, such as: maximum evaporation rates (1.5-4.4 mm/d), extinction depth (0.90 m), mineral dissolution (5.0e-9 g/d) and evaporation concentration (0.35 g/L). Despite the drawdown was uniformly distributed, the increase of groundwater salinity was rather uneven, while the temperature increase mimicked the atmospheric temperature increase. The initial groundwater salinity and the small changes in the evaporation rate controlled the evapoconcentration process in groundwater, while the effective porosity was the most sensitive parameter. This study demonstrates that shallow groundwater evaporation from sandy soils can produce homogeneous water table drawdown but appreciable differences in the distribution of groundwater salinity.


Introduction
The evaporation from bare soil driven by the presence of shallow groundwater is a significant component of the hydrological balance in semi-arid and arid locations (Assouline et al. 2014;Quinn et al. 2018;Shokri et al. 2010), although it could be relevant also in temperate climate, particularly in the Mediterranean area (Balugani et al. 2017). Direct evaporation from open waters have been studied in detail (Harwell 2012), as well evapotranspiration from reference crops (Jensen and Allen 2016), but less studies have been developed so far on bare soil evaporation induced by shallow groundwater (Allen et al. 2005;Bittelli et al. 2008;Flammini et al. 2018). The evaporation rate is primarily influenced by the microclimatic condition present at a given site, like temperature, wind speed, solar radiation, relative humidity, etc... (Allen et al. 2005;Martens et al. 2017); and by the soil's physical properties, like texture (Lehmann et al. 2018), soil organic matter content and water table depth (Alkhaier et al. 2012;Assouline et al. 2014). Coarse texture soils are characterized by high capacity to deliver water to the evaporation surface, but limited capillary raise due to lack of silt and clay fractions (Quinn et al. 2018). In fact, the liquid flux is driven by the capillary pressure gradients that develops in the vadose zone and remains active until the capillary gradients are larger than the gravitational and viscous forces (Lehmann et al. 2008). Thus, in coarse texture soils the evaporation rate is largely affected by the water table depth. An increase of the evaporation rate is usually followed by a soil temperature decrease because of the energy loss as latent heat flux (Todd et al. 2000). Many studies on the interaction between liquid water movement, water vapor transfer, and heat flux have found that the movement of heat and soil moisture are coupled (Kurylyk et al. 2019). Although thermal gradients affect the redistribution of water in soils, the most important process, which determines the coupling between water and heat, is the transport of latent heat by vapor flux within the soil and at the interface between the soil and the atmosphere. The soil water and energy fluxes from the land surface have been investigated in detail (Hingerl et al. 2016;Jin et al. 2019;Larsen et al. 2016;Trevisan et al. 2020), but very few studies considered the influence of shallow groundwater on soil surface temperature and the entire vadose zone (Alkhaier et al. 2009;Alkhaier et al. 2012;Doble and Crosbie 2017;Kollet and Maxwell 2008). Although valuable findings have been captured from these studies, most of them focused on the relationship between groundwater and soil surface evaporation. This shows that the influence of shallow groundwater on the evaporation process and its related impact on the subsurface energy balance along a soil profile is still an active area of research as pointed out by Trautz et al. (2018). The latter highlighted the importance of intermediate-scale experimentation (1-10 m) against the use of column scale data to derive generalizations about bare soil evaporation dynamic upscaling. Therefore, this study was conducted using a large tank with representative scales and instrumentations employed in the field to quantify the spatial distribution of evaporation from shallow groundwater in sandy soils. As the authors are aware, this is the first study that elucidates via numerical modelling the spatial distribution of the evaporation rate from shallow groundwater in well controlled laboratory conditions. This study also investigates the influence of surface evaporation on evapoconcentration effect (increased solute concentration due to evaporation) in both the vadose and saturated zones.

Experimental Set Up
A large tank (1.4 m x 4.0 m x 1.3 m), assembled with an internal structure of armed PVC and fastened on an external structure of natural wood enforced by steel scaffolding pipes, was filled with coarse sand materials (≅9.0 m 3 of sandy sediments and ≅0.5 m 3 of gravel). The sediments were excavated from a sand pit located along a meander of the Aspio river alluvial plain (Ancona, Italy). The relatively homogeneous nature of the sedimentary succession consists of coarsening-upward sandy sediments. The sediments were transported at the Hydrogeological Laboratory of SIMAU Department (Università Politecnica delle Marche), poured into the tank starting from the inflow gravel wall toward the outflow wall and compacted by a large shovel. The tank was filled to a height of 1.1±0.02 m. The natural compaction of sediment under saturated conditions was monitored for two months and after an initial localized collapse of ≅2.0 cm near the infiltration point the compaction was found to be negligible and the collapsed part refilled.
A constant head is present and can be applied to the tank via an external reservoir to create a steady state flux, but for this research was not used. The two gravel walls at the inflow and at the outflow were used to stabilize the groundwater flux (Giambastiani et al. 2013). The monitoring started in February 2020 and lasted for 4 months.
Eight piezometers with a bottom screen of 5 cm length and four fully screened wells were equipped with high resolution multi-level samplers (MLSs) and installed to form a semi-regular grid (Fig. 1). The wells were located along the central flowline of the tank and had an internal diameter (i.d.) of 5.0 cm, while the piezometers located between the wells and the side walls of the tank had an i.d. of 2.0 cm. Each MLS consisted of 6 HDPE tubes (4 mm i.d.) placed around the wells and piezometers, each HDPE tube was connected to a micro-screen of 0.5 cm length; the sampling ports were equally spaced every 10 cm, from 5 to 65 cm from the tank bottom. The MLSs were sampled for salinity measurement only once at day 110 to minimize the perturbation of the ongoing experiment, collecting 20 mL in each saturated port in the central transect C1, C2 and C3 (Fig. 1). The groundwater heads, temperature and electrical conductivity were monitored every 10 minutes using a Soil & Water Diver ® water level data-logger (Eijkelkamp, Giesbeek, The Netherlands), corrected for atmospheric pressure changes via a Barologger ® (Eijkelkamp, Giesbeek, The Netherlands) placed at the soil surface. The water level data resolution was 0.02% of the full scale and the accuracy was ±0.05% of the full scale. Three 5TE ® Meter probes (Meter Environment, Pullman, WA, USA) were installed inside the tank at 5, 20 and 40 cm below ground level (b.g.l.) to monitor volumetric water content (VWC), Temperature (T) and Soil Bulk Electrical Conductivity (ECb). Soil's ECb was converted in EC according to the model of Hilhorst (2000) and subsequently EC data were converted into salinity with standard conversion factors (APHA 2017). All probes were connected to a Meter data logger (ECH 2 O) recording every 10 min.
Physical parameters (grain size, bulk density, porosity, etc...) were calculated via dry sieving and gravimetric measurements on 5 randomly collected samples see Supplementary Information (SI) Table S1. The saturated hydraulic conductivity (k) distribution was estimated by slug tests (Bouwer and Rice 1976), via the Kozeny-Carman formula (Carrier 2003) and via a constant rate pumping test. The direct estimation of evaporation from the water table was calculated using the White (1932), detailed information are reported in SI (Estimation of evaporation from groundwater).  (Zheng and Wang 1999) for simultaneous solution of flow and heat transport equations so that the effects of variable density can be considered.
To simulate the groundwater drawdown and the heat pulses created by the changes in atmospheric temperatures, a tridimensional transient state fully saturated SEAWAT model consisting of 80 rows, 28 columns and 22 layers was used to get a uniform cell size of 5 cm 3 .
Before to start the experiment, the constant rate pumping test (4.0 l/min for 25 minutes) was analyzed via MODFLOW-2005 and PEST (Doherty 2010) by running a simulation with the Well package as sole boundary condition and using the Wetting factor to allow dry cells to be re-wetted once the pumping ceased. The automated inverse model PEST was used to obtain best estimates for the horizontal and vertical k values of the sand, the gravel walls and the average S y of the whole tank. The vertical k values were tied to the horizontal ones to constrain parameters estimation; the objective function consisted of piezometric heads collected every minute, recorded from 11 observation piezometers and the pumping well (B2 in Fig. 1).
All the boundaries were set as no flow boundaries except for the upper most saturated layer that was set as the evaporation surface by using the Evapotranspiration package of SEAWAT. The latter is defined by a maximum evaporation rate and an extinction depth from the ground surface.
To simulate the dissolution of soluble salts and mineral phases like carbonates, the mass loading rate package of SEAWAT was employed. The water table was uniformly set at 0.50 m from the ground surface before the start of the experiment, while the small differences in groundwater salinity and temperature were interpolated and used as initial conditions. The time was discretized into 29 stress periods of 5 days each, this allowed to mimic the temperature variation recorded by the 5TE probe at 40 cm depth by setting the mean values over 5 days in the Time Variant Specified Concentration package of SEAWAT. Flow and transport parameters used for the simulations are listed reported in SI Table S2.
The Geometric Multigrid (GMG) was used as numerical solver for the groundwater flow and the third-order scheme Total Variations Diminishing ULTIMATE (TVD) was used to solve the advective term of the solute transport problem. Courant number was set to 0.1 to further minimize numerical dispersion.
The results obtained from the calibration procedure were evaluated computing different statistical indices. In the present study, the robustness of the applied methodology was defined by means of three indices: coefficient of determination (R 2 ), Nash-Sutcliffe efficiency (NSE) and percent of bias (PBIAS). According to Moriasi et al. (2007), the satisfactory thresholds of the three statistical indices for an acceptable streamflow simulation are R 2 ≥ 0.50, NSE ≥ 0.50 and PBIAS ± 25%.

Unsaturated Zone Monitoring
In the unsaturated zone, the VWC showed a slow decline at all the monitored depths until day 115, where the topsoil started to desiccate more rapidly respect to the previous period (Fig. 2). The slow decline and not an abrupt one was imputable to the capillary rise forces, that provided a constant supply of water towards the evaporation surface, even if the sediment was a medium coarse sand. Congruently, the porewater salinity in all the monitored ports increased due to the evapoconcentration process (Mastrocicco et al. 2019;Tanji 2002) although the salinity did not increase linearly, but after a rapid rise it stabilized after day 60 and increased abruptly after day 120 in the topsoil.
This behaviour agreed with the monitored soil temperature, which oscillated around 20±2 °C until day 60 then slowly increased until day 120 and rapidly until the end of  Figure 3 shows the results of the inverse modelling of the constant rate pumping test, here it can be seen that the model was able to fully capture the heads variation in Fig. 3 Simulated hydraulic heads (red lines) and observed (black circles) in all the 12 monitoring points (1 per panel with their name in bold, refer to Fig. 1 for position) during the constant rate pumping test all the 12 monitoring points with a very low discrepancy between the observed and calculated values, exemplified by a R 2 of 0.981, a NSE of 0.98 and a PBIAS of only 0.08%. Apart from small overestimations and underestimations of the hydraulic heads in A2, A3 and C3 the model was able to accurately reproduce the heads drawdown and recovery. Although, it can be noticed that D1 observations are more scattered than the other ones; this was due to the pressure transducer that was factory calibrated for heads variations from 1 to 100 m (1 cm resolution), while all the others had pressure transducers calibrated for heads variations from 1 to 10 m, allowing a higher resolution (0.2 cm). Table 1 shows the observed mean k values and the respective standard deviations for the slug tests, here a slightly lower value was retrieved and a relatively high standard deviation. This is normal since the slug tests provide information in the vicinity of the monitoring point (Giambastiani et al. 2013;Paradis et al. 2016) and usually the values are lower than the k values provided by pumping tests (Butler and Healey 1998). This could be due to small scale variations of k or to partial development of the monitoring wells, although in this case it was most likely due to the first reason since the wells were fully developed. Moreover, the Kozeni-Carman equation produced a very similar average k value with an even higher standard deviation, witnessing that small k variations are pronounced. Table 1 shows that the horizontal and vertical k values for both the gravel wall and the coarse medium sand in the tank are within the range of typical values for such textures (Fetter 2001). Also the S y is within typical values and was unified for both sand and gravel since initial sensitivity analyses provided little changes using distinct values. Nevertheless, the S y resulted the most sensitive parameter, followed by the horizontal k of the sand, while the horizontal k of the gravel wall and the specific storage (S s ) were negligible. Vertical k parameters of sand and gravel were not considered in the sensitivity analysis since they were tied to horizontal k values.

K Field Characterization Via Inverse Modelling
From the hydraulic characterization of the tank the flow fundamental parameters were identified and employed to simulate the evaporation experiment. Figure 4 shows the observed and simulated trends of groundwater level, temperature, and salinity. The numerical model was able to capture the temperature variations even if the stress periods were set with a duration of 5 days; decreasing the stress period duration would allow to simulate more accurately the temperature variations during the experiment. Although, the model performance was considered satisfactory even using such a time discretization (Table 2), with all the three performance indicators well above the minimal requirement defined by Moriasi et al. (2007). It must be noticed that even if the initial temperature distribution was non-homogeneous, the time variant imposed temperature at the water table was set constant for the whole tank and this generated a nearly uniform temperature distribution in the numerical model (Fig. 4). This was also evident for the calculated heads, that were driven by the evaporation rate (the only stress set to the groundwater flow component of the model). This, however, was not held constant for the whole period, since after a few trials it was impossible to calibrate the model given the change in slope of groundwater drawdown recorded from day 115 (Fig. 4). Thus, the evaporation rate was split into 3 stress periods, the first from day 0 to day 60, the second from 61 to 120 and the third one from 121 till the end of the simulation, as identified by the vadose zone monitoring. The calibrated maximum evaporation rates and the extinction depth are reported in Table 2 (lower part). The rates are consistent with literature data on evaporation from bare soils in temperate climates (Aydin et al. 2008;Bittelli et al. 2008). Although the evaporation rates estimated via the White's method (see SI Fig. S1) were much lower than the surface evaporation rates and found to be similar to field observed values recorded in semi-arid climates (Gong et al. 2020).

Saturated Zone Monitoring and Modelling Performance
Surprisingly, the extinction depth found in this study, 0.90 m b.g.l. is higher than usually found in literature. In fact, Mansell and Hussey (2005) showed that when the water table is at 0.60 m b.g.l., evaporation from coarse and medium sands was approximately 10% of the evaporation from surface waters. The same conclusions were reached by Neal (2012). A review of the extinction depth values in response of soil texture and land cover was proposed by Shah et al. (2007) where for sandy soils they proposed the value of 0.50±0.05 m.
The model performance was very satisfying (Table 2 upper part) for groundwater level and temperature but the same was not true for salinity, which displayed lower values of both R 2 and NSE. This was likely due to the long screens of the monitoring wells that promoted artificial mixing (McMillan et al. 2014).
Besides, the evapoconcentration alone was not sufficient to describe and quantify the observed salinity increase, in fact without the minerals and salts dissolution, here accounted for with the mass loading rate, rather poor model performance indicators were obtained for salinity (not shown). The best model performance was obtained via trial and error and the values are reported in the lower part of Table 2. Figure 5 reported the concentration profile of calculated salinity at day 110 when the central MLSs (C1, C2 and C3) were sampled. An appreciable salinity gradient developed in all the MLSs although the groundwater was still fresh. Nevertheless, this proves that the increase in salinity observed with the continuous monitoring was essentially due to evapoconcentration processes and not due to salt and mineral phases dissolution. Moreover, this well explain the sudden salinity variations observed in the fully screened monitoring wells that could have been caused by small convective cells that often develop in open tube wells, induced by temperature gradients (Colombani et al. 2016;McHugh et al. 2012).
Given that the transport parameters are the most uncertain in this model (Table 2 lower part), a single parameter sensitivity analysis has been performed on the most critical ones, namely: dispersivity, mass loading rate and effective porosity.

Sensitivity Analysis
The dispersivity and mass loading rate were doubled and halved, while the effective porosity was increased or decreased of 10% to be consistent with the physically based values of this parameter (Fetter 2001). The results are shown in Fig. 6, here it is evident that the dispersivity did not affect appreciably the calculated salinity. While the mass loading rate, which resembles the minerals dissolution rate, moderately effected the calculated salinity especially in the scenario with 100% increase respect to the calibrated value, producing a general increase of approximately 50 mg/l at the end of the simulation. Finally, the effective porosity resulted to be the most influencing parameter, with calculated salinities that in some monitoring wells were nearly doubled in the scenario with the highest plausible value (0.22); or calculated salinities that did not vary during the simulation in the scenario with the lowest plausible value (0.18).

Conclusions
This study shows that evaporation from coarse medium sands can be relevant in temperate environments if the water table is near to the ground surface (0.5-1.0 m), with a measured extinction depth (0.9 m) higher than the ones proposed by previous Fig. 6 Salinity in the 4 monitoring wells and 8 piezometers for the calibrated model (blue lines) and in both the increased (red lines) and reduced (black lines) model scenarios studies with sandy sediments. This has also implications for laboratory experiments, like column or tank experiments, where this important parameter is often unaccounted. Moreover, the monitoring data showed that groundwater temperatures and levels were homogeneously distributed over the monitoring network and their behavior could be approximated using a one-dimensional approach, treating the whole tank as a homogeneous soil column. Despite the relatively homogeneous drawdown, the White's method used to calculate the daily evaporation rate from the water table produced values affected by large spatial variability. This, in conjunction with small differences in the initial concentrations, determined appreciable changes in the evapoconcentration process that were well captured only using a three-dimensional flow and transport model. The sensitivity analysis showed that small changes in the effective porosity could lead to very different concentrations' distribution during the evapoconcentration process, thus the influence of this parameter should be carefully characterized in future studies. This study demonstrates that shallow groundwater evaporation from sandy soils can produce homogeneous water table drawdown but appreciable differences in groundwater salinity at the meso-scale and to accurately capture such variations high resolution multi-level samplers must be employed.