Sensitivity of Vadose Zone Water Fluxes to Climate Shifts in Arid Settings Open Access

In most arid countries groundwater is the primary freshwater source since surface water resources are limited and temporally extremely variable, resulting from highly temporal fluctuations in precipitation patterns at annual scale (Littmann and Berkowicz, 2008; Scanlon et al., 2006). With increasing population and industry causing an increase in water demand, groundwater resources are prone to depletion, and water scarcity is a severe issue for arid regions (FAO, 2003). To evaluate groundwater availability in the future, one major task is to estimate groundwater recharge. Strategies to estimate groundwater recharge in arid regions include isotope fractionation due to seasonal variability in precipitation and evaporation, water balance modeling, chloride-mass balancing (Subyani and Şen, 2006), the quantification of infiltration by lysimetry, tracer tests, and empirical approaches (see overviews given in, e.g., Gee and Hillel, 1988; de Vries and Simmers, 2002; Flint et al., 2002). Generally, groundwater recharge by direct infiltration and percolation is assumed to be low in the case of high evaporative losses. Nevertheless, especially in regions covered with highly permeable dune sands, precipitation can rapidly percolate to depths not affected by evaporation, and contribute to groundwater recharge (Dincer et al., 1974; Scanlon et al., 2006). As the unsaturated dune sands vertically can extend up to hundreds of meters, vadose zone processes governing the up- and downward movement of the infiltrating water have to be considered for estimating direct recharge rates in these areas (Scanlon et al., 1997). Efforts have been made to better understand and quantify vadose zone processes and water fluxes by improved measurement techniques of water content (e.g., Dahan et al., 2003; Jones et al., 2005), observation of scale effects (e.g., Hendrickx and Flury, 2001; Hopmans and Schoups, 2006; Mattson et al., 2004), and numerical modeling (e.g., Dong et al., 2003; Saito et al., 2006; Sakai et al., 2009). However, especially for arid regions, there are still many uncertainties due to difficult field conditions, such as highly variable or sporadic precipitation events, extreme evaporation rates, and high temperature variations between day and night, as well as coupled processes of water, vapor, and heat transport (Heitman et al., 2008).

In numerical investigations, the influence of temperature gradients as well as ambient soil water contents close to residual saturation have to be considered carefully if nonisothermal water fluxes are simulated in the vadose zone under arid conditions. Additionally, accurate calibration data obtained under nonisothermal conditions, including transient soil moisture variations, may help to reduce model uncertainties (Inoue et al., 2000). However, soil hydraulic characteristics governing unsaturated water flow are often measured under multistep equilibrium and isothermal conditions (Wildenschild et al., 2001) and are not appropriate under nonisothermal conditions at field scale. Inverse modeling strategies, including sensitivity and uncertainty analyses, and parameter estimation procedures of large-scale soil column experiments under controlled boundary conditions can notably improve the reliability of model predictions (Finsterle, 2004).

Besides the hydraulic soil parameters, ambient climatic conditions are relevant to analyze infiltration and percolation processes in the vadose zone of arid regions. Due to improved and enhanced climate modeling strategies, recent research studies have focused especially on scarce groundwater systems that result from decreased infiltration, which is driven by changes in temperature and precipitation during the past 100 years (e.g., Bou-Zeid and El-Fadel, 2002; Kundzewicz and Döll, 2009; Taylor et al., 2012). Climate change modeling concepts are based on greenhouse gas emission scenarios as proposed by the Intergovernmental Panel on Climate Change (IPCC, 2001). For many semiarid and arid regions, and specifically for North Africa and large parts of the Arabian Peninsula, increased temperatures and decreasing precipitation amounts during the summer months are expected for 2100 (Ragab and Prudhomme, 2002). The dependency of groundwater recharge on these climate projections can provide valuable information about long-term changes in the renewability of freshwater resources in regions that are already under severe water stress. Studying the impact of climate changes on groundwater recharge was recently the focus of large-scale hydrological and agricultural studies (e.g., Goderniaux et al., 2011; van Roosmalen et al., 2009). A specific focus of groundwater resources under climate change pressures in arid regions was set by Engelhardt et al. (2013a,b) and Engelhardt and Prömmel (2012). Nevertheless, such large-scale studies do not investigate the processes specific to the vadose zone. Reasons for neglecting vadose zone processes in large-scale groundwater recharge models are mainly increased computational requirements and the difficulty of inferring the soil hydraulic parameters for unsaturated zone flow calculations and the subsequent parameter estimation requirements (Vrugt et al., 2004). Currently, investigations coupling long-term regional shifts in climate with process-orientated, small-scale vadose zone models are still lacking (Candela et al., 2012). This may be attributed to the challenge of downscaling the results of climate model projections suitable for simulations of up- and downward water and water vapor fluxes at the more local scale of the vadose zone. Thus, not only climate and regional hydrological patterns have to be considered to study recharge in arid regions properly, but also, water, vapor, and heat fluxes in the context of hydraulic and thermal soil properties must be specified.

In this study, we present an inverse model of infiltration experiments accounting for variably saturated nonisothermal water fluxes and subsequent scenario modeling linking the results of the inverse model with projections of a global climate model. A study area in the eastern and currently arid to hyper-arid part of Saudi Arabia was selected. The goals of our research were (i) to identify physical soil parameters of primary importance regarding groundwater recharge in arid regions, (ii) to determine the most important atmospheric parameters for direct groundwater recharge, and (iii) to quantify the effect of different climate change projections based on greenhouse gas (GHG) emission scenarios on groundwater recharge in nonvegetated sand dune areas in an arid environment. The inverse model estimates hydraulic and thermal parameters and assesses their sensitivity for simulating heat, water, and vapor fluxes under transient arid climate conditions. The scenario modeling couples global climate model projections with regional recharge estimates. This allows insight into the influence of changing climatic conditions, as well as to derive threshold values for infiltration, evaporation, and percolation rates. Three GHG emission scenarios (A2, A1B, and B1) (Nakićenović and Swart, 2000) that give global climate projections for 2100 were investigated. They were simulated with the coupled atmosphere–ocean model EGMAM (Huebener et al., 2007). This allows transferring the influence of expected greenhouse gas emissions into future changes in groundwater recharge. Thus, a holistic analysis of hydraulic processes in deep vadose zones in response to shifts in climate in a present-day arid setting is obtained.

The investigated arid setting is located in the eastern part of the Arabian Peninsula and is bounded by the latitudes 29°75′ and 22°25′ N and longitudes 46°88′ and 54°38′ E (Fig. 1). The boundaries of the study area were chosen (i) to capture one of the biggest aquifer systems of the world, the Upper Mega Aquifer System, (Margat, 2007), (ii) to investigate an arid environment with thick vadose zones above the aquifers, and (iii) to extend the research area over several grid cells of the climate model.

The climate of the investigated setting is mainly arid with low precipitation and high inter-annual variability in rainfall (NOAA, 2010). Generally, precipitation only occurs during the winter months and dry and hot summers prevail. In the investigated region, the long-term average annual precipitation measured from 1978 to 2012 ranges between 90 110 mm yr⁻¹, whereas in the southern parts annual precipitation can decrease to zero. In the northern part of the study site in particular, single precipitation events can amount up to 200 mm (Engelhardt et al., 2013b). Nevertheless, on an annual basis, potential evapotranspiration by far exceeds precipitation in the whole region, with values up to 4000 mm yr⁻¹ due to the dry and hot summer climate (ElNesr et al., 2010). In the study area, average annual temperatures vary about 25°C, with maximum daily values of up to 47°C during summer, dropping down to 3°C in winter (Engelhardt et al., 2013a,b).

Due to the low precipitation, large parts in eastern Saudi Arabia show a reduced soil development without a clear soil structure within the profile (Shadfan et al., 1984). In large parts of the region above the aquifer thick unconsolidated Quaternary sediments are deposited, forming extensive sand deserts and gravel sheets (Shadfan et al., 1984). The large dune regions offer a high potential for infiltration and direct groundwater recharge (Dincer et al., 1974, Engelhardt et al., 2013b).

For the numerical investigations, small column experiments of hydraulic properties from repacked soil cores were used to determine initial estimates of hydraulic parameters for inverse modeling studies. The parameter calibration and a sensitivity analysis were performed with drainage and evaporation data collected from large column experiments with 1-m length. Additionally, long-term predictive modeling was conducted to assess the effect of different GHG emission scenarios and corresponding predicted changes in climate on downward water fluxes in the deep vadose zone of the investigation area in south-eastern Saudi Arabia.

Large column experiments (1-m length) were conducted with repacked homogenous red dune fine sand (100% sand content) that was sampled near Riyadh (25°13′ N and 47°61′ E). Hydraulic and thermal properties of the dune sand were determined before the large column experiments in small column experiments at the Environmental Molecular Science Laboratory (EMSL, Richland, WA, USA). These measurements were performed on repacked 100-cm³ samples. The samples were packed to a bulk density of 1700 kg m⁻³ as measured in the field. Saturated hydraulic conductivity was measured with a constant head permeability test (Wietsma et al., 2009). The experimental water retention curve was obtained by a hanging column apparatus with a 100-mL probe and 15 measurement levels (10, 20, 30, 40, 45, 50, 60, 70, 80, 90, 100, 125, 150, 175, and 200 cm). The measured water retention data was used to estimate the soil hydraulic parameters of the Brooks and Corey (1964) model with the program SWRC Fit (Seki, 2007) (Table 1). The Brooks and Corey model was chosen due to its validity for coarse grained material (e.g., Schaap et al., 2001). The derived parameters were used as initial parameter estimates for inverse modeling using discharge and evaporation data from large column experiments. Retention measurements were performed up to a suction of −200 cm. Fitting the Brooks and Corey model to the measurements resulted in an estimated residual water content of 0.067 cm³ cm⁻³. Nevertheless, residual water content can be better defined as remaining water content at high capillary pressures. Thus, another fitting of the Brooks and Corey model to the retention data was performed with the residual water content restricted to zero and the capillary pressure to 10^4.2 cm, resembling wilting point (Groenevelt and Grant, 2004). High capillary pressures are expected to occur under arid conditions, and setting the residual water content to zero was expected to reflect particularly dry conditions better. This was also expected to provide more suitable initial hydraulic parameters for the inverse modeling.

Thermal conductivity was measured by a transient line heat source method with a KD2 thermal analyzer (Decagon Devices) for five steps from dry to fully saturated conditions. The 30-mm needle sensors were installed in the 100-cm³ repacked samples. The sample was homogenously desaturated for each measurement step. Nevertheless, measurements of thermal conductivity are linked to several uncertainties, especially when relating them to low saturated conditions under variable temperatures. Uncertainties are due to differences in sample temperature, influences of the measurement device, and nonequilibrium conditions caused by diffusive effects (Smits et al., 2013).

For the large column experiments, the dune sand was packed with a bulk density of 1700 kg m⁻³ and a water content of 0.07 m³ m⁻³ was adjusted. In the column with a diameter of 0.19 m and a height of 0.98 m, including a 0.05-m controlled air head space at the top and a bottom suction plate, measurements were made for changes in water content and temperature. The air-filled head space at the top of the column was used to generate a flow of air, which was controlled in velocity and relative humidity applying a flowmeter and forcing the air flow along a specific saturated salt solution to constrain the humidity. By measuring the relative humidity and the temperature of the inflowing and the outflowing air, evaporation and soil water uptake were calculated. The suction plate enabled column outflow with adjustable pressure head conditions that were set lower than atmospheric pressure to prevent water accumulation at the column bottom. Discharge passing the suction plate was measured continuously. Irrigation, which was controlled in intensity and total amount, was applied at the column top. Additionally, capillary spirals, which were heated or cooled by circulated water, controlled the temperature both on top and at the bottom of the column. Water content and temperature profiles, as well as water fluxes at the column boundaries, yielding evaporation and discharge, were measured at high temporal resolution of up to 1 min.

The inverse modeling of soil hydraulic and thermal properties using large column experiments was based on the measurements of transient water contents and temperatures at depths of 0.04, 0.14, 0.24, 0.34, 0.44, 0.60, 0.70, and 0.84 m, as well as transient flux measurements of evaporation at the top of the column and discharge at the bottom with a temporal resolution between 1 and 30 min, respectively. Total duration of the experiment was 300 h. The experiment included two irrigation events: an initial period of 5 h with a total of 170 mm until discharge was obtained followed by 120 h without irrigation and a second irrigation period of 2 h with 25 mm. These irrigation patterns were chosen (i) to have an initial percolation throughout the whole column for analyzing wetting processes and (ii) for analyzing percolation and infiltration depths of lower irrigation amounts. The temperature at the top of the column was set to 40°C during the first 200 h to simulate hot summer events and decreased to 35°C for the remaining 100 h to evaluate small temperature changes for soil surface water fluxes. At the bottom of the column temperature was fixed to 18°C to mimic the neutral temperature zone. The laboratory column experiments are described in detail in Pfletschinger et al. (2012).

To take into account the changing climate in Saudi Arabia until the end of the century, global simulations of future climate were analyzed. The global climate simulations used in this study were performed within the EU-project ENSEMBLES (http://www.ensembles-eu.org) with the global climate model EGMAM (ECHO-G Middle Atmosphere Model; Huebener et al., 2007), which is a vertically extended version of the coupled atmosphere ocean model ECHO-G (Legutke and Voss, 1999). It consists of the atmosphere model ECHAM4 (Roeckner et al., 1996) and the ocean model HOPE-G (Wolff et al., 1997). In the EGMAM setup ECHAM4 supports a horizontal resolution of approximately 3.75° × 3.75°. The troposphere and stratosphere up to approximately 80 km are represented by 39 levels. HOPE-G has 21 layers and a horizontal resolution of approximately 2.8° × 2.8° with a grid refinement near the equator up to 0.5°. For coupling, a flux correction is applied for heat and freshwater exchange. The orography of the grid cells is based on a digital elevation model that is averaged for the grid cells of the climate model to a horizontal resolution of 3.75° × 3.75°.

Future climate strongly depends on past and future emission of greenhouse gases, which is used as forcing for the climate simulations. For the spin-up, from 1860 to 2000, a 1% yr⁻¹ CO₂ increase experiment, called dCO₂, is used. This simulation starts in 1860 with a CO₂ concentration of 280 ppm that increases each year by 1%. A doubling of the initial CO₂ concentration is reached in 1930. Afterward the concentration is kept constant until 2000.

Different scenarios of future greenhouse gas emissions were developed for the IPCC Third Assessment Report, the IPCC Fourth Assessment Report, and described in Nakićenović and Swart (2000). The scenarios are based on various storylines describing the possible developments of economy, population growth, energy use, technology development, land-use, and policy strategies. For this study, three scenarios were chosen that represent the possible range of future greenhouse gas emissions, but without being borderline scenarios. Scenario B1 describes a convergent world with emphasis on global solutions to economic, social, and environmental sustainability (Nakićenović and Swart, 2000). Scenario A1B describes a world with economic and cultural convergence and capacity building, with a substantial reduction in regional differences in per capita income, and the rapid introduction of new and more efficient technologies with a balance across fossil and non-fossil energy sources (Nakićenović and Swart, 2000). Scenario A2 describes a heterogeneous world with a strengthening of regional cultural identities, with an emphasis on family values and local traditions, high population growth, and less concern for rapid economic development (Nakićenović and Swart, 2000).

For each GHG scenario three realizations were simulated that describe similarly probable conditions (Huebener et al., 2007). We randomly selected one realization from each scenario to preserve the day-to-day variability. The use of the ensemble mean over the three realizations would have resulted in an underrepresentation of the climate variability and therefore in a loss of extremes such as intense precipitation events or dry periods. To convert the climate simulation results into input data for the simulation of vadose zone fluxes, grid boxes covering the study area have to be extracted. The corresponding three selected grid boxes (48.75°E–27.83°N, 48.75°E–24.12°N, 52.5°E–24.12°N, Fig. 1) represent the low-lying arid to hyper-arid areas in eastern Saudi Arabia. For the vadose zone scenario modeling, daily values of precipitation, surface temperature, minimum and maximum temperature, relative humidity, wind speed 10 m above the surface, and surface solar radiation are averaged over the three grid boxes and assigned as upper atmospheric boundary conditions.

The governing equations for heat transport consider the conduction of sensible heat, convection of sensible heat by liquid water and water vapor, and the convection of latent heat by vapor flow (Šimůnek et al., 2009). The hydraulic and thermal equations are solved in a fully coupled way.

The inverse model was set up to estimate the hydraulic soil properties from column experiments (Fig. 2). The 0.93-m soil column was spatially discretized with 931 nodes at a constant distance of 0.001 m. Temporal discretization was set to 3.6 s.

The hydraulic boundary on top of the column accounted for hourly changes of precipitation and potential evaporation. The measured transient temperature function was assigned to the top of the column for the soil temperature in daily resolution. At the top, the system switches between a flux and a head boundary condition, constrained by the minimum head at the uppermost node, until which potential evaporation equals actual evaporation. As soon as this maximum head, which was set to −2 m, is exceeded, the actual evaporation rate decreases and is constrained by the availability of water in the soil (Šimůnek et al., 2009).

At the bottom, the temperature was fixed to 18°C. A suction plate was installed at the bottom of the column, which was represented in the model by defining hydraulic properties different from the soil material for a 0.1-m-thick bottom layer. Hydraulic properties of the suction plate were initially defined as given in the manufacturer’s data sheet (Soil Moisture Equipment Corp., 2009), but later adjusted in a joint inversion with the hydraulic parameters of the dune sand. Thermal properties of the suction plate were set equal to the dune sand as the heating element was placed above the suction plate. In the experiments, the pressure head was set to −0.6 m at the suction plate to enable downward flow. A seepage face was assigned to the column bottom to mimic the suction plate. The seepage face boundary implies that zero-flux boundary conditions prevail as long as the local pressure head at the bottom of the profile is below the fixed seepage face pressure head value. As the local bottom pressure head reaches the defined seepage face value, the boundary discharge is constrained by the prescribed pressure head. Due to uncertainties in maintaining the pressure head exactly at −0.6 m at the bottom of the column during the experiments, this value was included in the parameter estimation process of the inverse model.

Initial conditions were specified in terms of the measured water content and temperature distribution at the beginning of the simulation due to the experimental setup with a homogenous water content of 0.07 m³ m⁻³ and a homogenous temperature of 20°C. In the experimental setup homogenous initial conditions were ensured by stepwise packing and water content and temperature monitoring during and after packing (Pfletschinger et al., 2012).

Automated model calibration was performed with the nonlinear parameter estimator PEST (Doherty, 2010) through an inverse procedure based on the Gauss–Marquardt–Levenberg method. The parameter optimization process in PEST uses model outputs and corresponding observation measurements by minimizing the weighted sum of squared differences. Weighting factors were assigned to increase the impact of the water and temperature breakthrough on the model calibration and to capture possible measurement errors. Weighting factors of breakthrough measurements were set to 10. Generally, weighting factors for temperature measurements were reduced to 0.1, which accounts for the higher numbers of temperature values compared to the small numbers of water content and fluxes. PEST optimization includes computing sensitivities of all observation points, providing a measure of how much a simulated value changes in response to the perturbation of an adjustable parameter. Composite sensitivities are obtained for all observation data and estimated model parameters. In the automated calibration 762 observation data were included, given by measured discharge at the bottom of the column; evaporation at the top of the column; temperature continuously measured at 0.04, 0.14, 0.24, 0.34, 0.44, 0.60, 0.70, and 0.84 m with a temporal resolution between 1 min and 10 h; and volumetric water content continuously measured at 0.04, 0.14, 0.24, 0.34, 0.44, 0.60, 0.70, and 0.84 m with at temporal resolution between 1 min and 20 h.

With the inverse model the following hydraulic model parameters of the dune sand and the suction plate were estimated: the residual volumetric water content θ_r, the saturated volumetric water content θ_s, the Brooks and Corey parameters h_b and λ (Eq. [2]), and the saturated hydraulic conductivity K_s (Eq. [3]).

Due to their low sensitivity and to avoid an overparametrization the tortuosity l of the hydraulic properties of the dune sand and suction plate and the three empirical parameters describing thermal conductivity of the dune sand as function of volumetric water content were fixed to the value derived from the laboratory measurements of the soil retention curve and thermal conductivities, respectively.

For the scenario modeling, the model domain was extended to 30 m to resemble a vadose zone above the groundwater surface. Spatial discretization was changed to a resolution of 0.05 m, and the temporal discretization ranged between 0.05 and 1 d. At the bottom, a free drainage boundary condition was defined to account for a deep groundwater surface beneath the vadose zone. The bottom temperature was fixed at 18°C to mimic groundwater temperature.

A spin-up period was run before the future scenario modeling to create initial conditions that resemble the effect of long-term arid conditions on the soil water profile. For the spin-up period, monthly varying climate data obtained from the EGMAM simulations for the period 1860 until 2000 based on the dCO₂ experiment were assigned to the top of the model. The initial water content at the beginning of the spin-up period was set to 0.06 m³ m⁻³ to account for an initially wetted profile. Temperature distribution was set homogenously to 18°C.

For the scenario modeling of the vadose zone water fluxes from 2000 until 2100, boundary conditions at the top of the model were defined with daily variable input data. Atmospheric boundary conditions were extracted from the EGMAM simulations and captured precipitation and mean air temperature. The potential evaporation is calculated in Hydrus with the Penman–Monteith combination equation (Monteith, 1981; Šimůnek et al., 2009). For this, daily input data for minimum and maximum air temperature, relative air humidity, wind velocity, and monthly averaged solar radiation, also extracted from the EGMAM simulations, were assigned as climate conditions between 2000 and 2100.

Mean residuals for water content were 0.02 m³ m⁻³ (0.63 m depth), 0.025 m³ m⁻³ (0.23, 0.33, 0.43, 0.53, 0.73 m depth), 0.03 m³ m⁻³ (0.13 m depth), and 0.035 m³ m⁻³ (0.03 m depth). The higher residuals at the top of the column resulted from the first measurements after the beginning of the irrigations. Despite single high residuals between measured and calculated discharge at the beginning of the experiment, overall discharge volumes and evaporation amounts showed a good match indicated by mean residuals of 0.005 and 0.0003 m, respectively (Fig. 3). In the experiments, some fluctuations in the applied suction at the bottom, as well as in the incoming air relative humidity at the top of the column accounted for small variations within the measured evaporation that were not simulated, giving residuals for single evaporation measurements of up to 0.1 m. For the temperature measurements, maximum residuals were up to 8°C for single values at 0.04 and 0.14 m depth during initial heating. To avoid overparameterization of the inverse model, thermal parameters that were assessed by PEST with very low sensitivities were excluded from the final optimization runs.

For the dune sands, the highest dimensionless scaled sensitivity coefficients were calculated for air-entry pressure. The air-entry pressure defines the pressure head at which the water content drops below θ_s and hydraulic conductivity drops below K_s as air enters the pores that are formerly completely filled with water. This high sensitivity could probably be explained by the fact that the air-entry pressure shifts the entire retention curve and the unsaturated hydraulic conductivity function. Especially in sandy soils, water retention and hydraulic conductivity significantly decrease after air enters the porous structure, indicated by the steep slope of the retention curve (Fig. 4). For the fine dune sand used in our study a rather high air-entry pressure of 50 cm was measured in the small column experiments (Fig. 4) compared to investigations analyzing retention data of other dune sands (e.g., Toride et al., 2003; Sakai et al., 2009). The inverse modeling confirmed the range of the air-entry pressure for the large column experiments with an optimized value of 31.5 cm. For the large column experiment, the air-entry pressure is of high importance, as it defines the change from fast downward water drainage to capillary driven decreasing water fluxes. Compared to small-scale retention measurements, in the large column experiments the air-entry pressure controls the water movement in the whole column given by up- and downward fluxes that are influenced by this decrease.

The residual and saturated water contents were assessed as medium sensitive. Low sensitivities were calculated for K_s and the pore size distribution coefficient λ of the Brooks–Corey soil model. The low sensitivity of both parameters probably reflects the rather dry conditions that prevailed during the experiments whereas most of the studies which detected a high sensitivity of K_s and λ rather focus on the wet part of the retention curve. For the suction plate, dimensionless scaled sensitivity coefficients were at least one third lower compared to the dune sand parameters. Highest sensitivities were also calculated for the Brooks–Corey parameter h_b, while the other parameters obtained sensitivity coefficients near zero.

Overall, the retention curve obtained from the inverse model of the laboratory experiments shows a shift along the water content line due to an increased saturated water content (Fig. 4). At volumetric water contents below 0.08 (m³ m⁻³) the pressure head measured with the small columns rises to infinity and water flow is disabled, while the pressure head derived from the parameter estimation with PEST increases at low volumetric water contents more smoothly, and water flow is still possible. Thus, distinct differences between the retention curve derived from the small column and large column experiments are obtained close to the residual water content. These differences mainly result from the different imposed boundary conditions. The small column experiments were performed under stationary and isothermal conditions, whereas the large columns were performed under nonisothermal and highly transient boundary conditions.

Regarding observation data used in the calibration process, sensitivities were highest for water content observations at the time of breakthrough of the infiltration water in the corresponding depths, whereas top observations had higher sensitivities with respect to deeper observations. During gravitational flow, sensitivities declined, and increased again when fluxes decreased. For these two stages, bottom observations had higher sensitivities than top observations. During the drying phase, where water fluxes converged to zero, sensitivities decreased again, whereas sensitivities at the top were higher presumably due to upward evaporative fluxes. Therefore, if monitoring campaigns are designed after precipitation events, shallow parts of infiltration fronts are most important to monitor and will give highly sensitive data for numerical investigations. However, in an arid environment, many days without precipitation prevail. During these periods, deeper parts of the soil water profile will also provide sensitive and valuable data sets to simulate vadose zone water fluxes in an arid setting.

The climate models compute an overall increase in temperature in the study area by approximately 2°C, with the highest temperature increases by 4 to 6°C for the A2 scenario and the lowest temperature increases of 1 to 2°C for scenario B1 (Fig. 5). Lowest total cumulative precipitation within the 100-yr scenario period was obtained for the selected A1B realization with 8942 mm, whereas highest total precipitation was given for the selected realization of scenario B1 with 10,108 mm. The total difference of the cumulative precipitation between the three realizations remains only minor at 1166 mm. The frequency distribution of daily precipitation is very similar in all scenarios, showing highest frequency of no precipitation at all, and quickly decaying frequencies with growing precipitation amounts. However, Scenario B1 tends to have fewer days without precipitation and fewer days with very high precipitation amounts, but more days with medium precipitation amounts compared to Scenarios A1B and A2. This leads in total to the slightly higher cumulative precipitation amount. In all scenarios the annual variation in precipitation is quite strong with highest annual sums of 160 to 220 mm yr⁻¹ (Fig. 5). This is also visible in the time-series of daily precipitation sums with highest single precipitation events in all scenarios of 40 to 55 mm d⁻¹. The analysis of the different realizations of each scenario demonstrates that the occurrence of these high precipitation events is randomly distributed over time without favoring a special period. Additionally, there is no significant trend in precipitation. Therefore, changes in the GHG concentration do not strongly force precipitation over Saudi Arabia, but mainly temperature.

The strong variation in daily as well as annual precipitation and the differences between the realizations for all three scenarios illustrate the high internal variability, and therefore the high uncertainty and noise, of simulated precipitation over Saudi Arabia. There is also no common trend identified comparable to that of temperature. This is valid for all widely used global climate models due to uncertainties in hydrological response (Christensen et al., 2007).

During the spin-up period, from 1860 to 2000, the water content in the 30-m-deep profile converged to steady-state conditions with a water content of 0.047 m³ m⁻³ which deceased to lower water contents down to 0.02 m³ m⁻³ in the upper 2 m. The water content below 2 m reflects field capacity conditions, where water is held against gravitational potential. In the upper 2 m, water content is additionally influenced by the atmospheric conditions. Here, depending on precipitation amounts and atmospheric pressures that influence evaporation, water fluxes are either directed down- or upward. At the end of the spin-up period, the discharge reached steady-state conditions with 13 mm yr⁻¹. With respect to groundwater recharge derived for the study area from hydrological simulations at catchment scale by Engelhardt et al. (2013b), resulting in about 5% of the precipitation acting as recharge, the obtained recharge value seems to be relatively high, showing that 12 to 18% of the annual average precipitation contributes to groundwater recharge. Limiting factors such as soil cover, spatial variation of precipitation, land types, and topography are disregarded in the numerical analysis of the soil column experiments and potentially decrease recharge volume.

For all three GHG emission scenarios, in 2100 the mean soil water storage was around 1500 mm in the whole 30-m soil profile and mean recharge rates varied around 22 to 29 mm yr⁻¹. Soil water storage and bottom flux, and thus recharge rates, are highly affected by precipitation pattern. The peaks in soil water storage are driving higher discharge rates and can be related to high annual precipitation values (Fig. 6) mainly resulting from high single precipitation events or a fast sequence of following precipitation events. The highest recharge peak was observed for the B1 scenario with 55 mm yr⁻¹ emerging from a high annual precipitation of more than 250 mm in the year 2037. Several periods with small recharge rates of 13 mm yr⁻¹ were also observed for the B1 scenario. For the A1B scenario, maximum and minimum recharge rates were 51 and 13 mm yr⁻¹ in the year 2045 and 2030, respectively. In the A2 scenario, maximum recharge was obtained in 2075 with 51 mm yr⁻¹, resulting from a high annual precipitation in 2071 of around 140 mm with one high precipitation period of 70 mm within 5 d (Fig. 7). Minimum recharge during the A2 scenario was calculated with 15 mm yr⁻¹ for 2025. In all scenarios, maximum recharge is driven by single precipitation events occurring between 2037 and 2075. Nevertheless, as the precipitation patterns of the climate scenarios are highly variable for different realizations, attributing future recharge patterns to time or years rather than to precipitation amounts, is quite vague. All scenarios also predict that future minimum recharge rates will only slightly decrease below current recharge rates. Cumulative recharge during the 100-yr simulation time is highly similar for all the three scenarios. Total cumulative recharges in the year 2100 are 2200 mm for A1B, 2310 mm for A2, and 2320 mm for B1, accounting to 25, 29, and 23% of the cumulative precipitation of the whole 100-yr simulation period.

In terms of cumulative recharge values, the soil water storage is rather similar within the whole 30-m-deep profile for all three scenarios. The absolute value of soil water storage gradually increases from around 1.4 m³ to around 1.55 m³ in the whole soil profile from 2000 to 2100. This is contributed by the daily precipitation values that account for lower evaporation losses due to deeper infiltration used in the predictive modeling from 2000 to 2100, compared to averaged monthly precipitation values assigned during the spin-up period. This emphasizes the need of climate prediction at high temporal resolution if groundwater recharge should be estimated as a response to future climate changes.

Water content changes in the whole soil profiles for all three scenarios basically show similar percolation patterns as a response to high precipitation events, even though the timing of the single events differs. Percolation fronts of a fast sequence of high precipitation events can be traced in all three scenarios down deeper parts of the vadose zone. However, probably due to lower precipitation intensities the infiltration fronts in A1B only reach down to a depth of 15 m, while in in B1 and A2 high precipitation events are retrieved down to a depth of 30 m (Fig. 7).

When zooming in to the upper 5 m of the soil profiles, an evaporation depth of less than 3 m can be observed for all three scenarios. Even the larger increase in temperature in the A2 scenario does not show significant increases in evaporation depth at the end of the simulation period (Fig. 8).

Comparing water and vapor fluxes in summer and winter show different flux patterns driven by temperature (Fig. 9). In summer time, in the deeper part of the vadose zone, both water and vapor fluxes are mainly directed downward. Only in the upper few centimeters of the soil profile, vapor flux is directed upward, contributing to evaporation. During winter, the vapor flux is directed upward down to 1 m soil depth, which is driven by thermal vapor fluxes, as the thermal gradient is negative and thus directed upward. Water fluxes during winter are directed downward in the whole soil profile, as precipitation mainly occurs in the winter time leading to gravitationally driven downward water flow. Below the first centimeter of the soil profile, both summer and winter water fluxes are directed downward. In the deeper part of the soil profile (blow 2 m depth) summer water fluxes approach zero, while winter water fluxes reach steady-state conditions. Thus, only during winter month infiltrated precipitation will reach deeper parts of the vadose zone in higher amounts and thus can enable groundwater recharge.

The inverse modeling based on sensitivity coefficients revealed that highly relevant parameters for groundwater studies in arid regions are air-entry pressure and residual water content, rather than saturated hydraulic conductivity. Thus, measurements of the retention curve near residual saturation can be more helpful for characterizing water flow patterns in arid regions than for conducting infiltration tests. Nevertheless, measuring retention curves in the small columns still might not reflect the retention curve valid under field conditions. Transferring the results obtained by soil column experiments to field scale still lacks in taking heterogeneities, as well as other environmental factors like vegetation, for example, into account. Thus, additional experiments employing large-scale lysimeter setups under real field conditions, including event-based soil moisture monitoring are highly recommended for future research.

The predictive modeling studies based on GHG emission scenarios showed an increase of the mean recharge by 7 to 25 mm yr⁻¹ in 2100. However, obtained values only account for homogenous high permeable sand profiles and effective recharge rates are expected to be under field conditions distinctively smaller. Regarding climate, highest influences for groundwater recharge are driven by precipitation oscillations that occur in all three climate predictions, A1B, A2, and B1. Even the high temperature increase in Scenario A2 does not affect groundwater recharge rates. For groundwater recharge predictions in arid regions, reliable precipitation predictions, especially on a daily basis, are most relevant to reduce uncertainty in numerical investigations.

Due to the uncertainty and noise of the precipitation in the climate simulations it is not possible to predict the particular time of high precipitation events. However, the simulations with the climate model indicate for all GHG scenarios that there will be periods in future with higher precipitation events leading to increasing soil water storage, probable groundwater recharge, and a refreshing of the groundwater resources.

The study showed that groundwater recharge is mainly governed by precipitation rates, in which single high precipitation events, which can deeply infiltrate, are the main source for direct recharge in an arid environment. Different climate scenario studies corresponding to different GHG emission scenarios showed only marginal effects on recharge rates, and thus cumulative groundwater recharge rates until 2100 vary only minor for different climate predictions.

We highly acknowledge the soil hydraulic and thermal property measurements of the dune sand performed at the Environmental Molecular Science Laboratory (EMSL, Pacific Northwest National Laboratory, Richland, WA, USA) with special thanks to M. Oostrom and T. Wietsma. The work was founded by the BMBF (German Federal Ministry of Education and Research) within the IPSWaT (International Postgraduate Studies on Water and Technologie) program. Climate simulations were performed and founded within the EU project ENSEMBLES (http://ensembles-eu.metoffice.com/).