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

In arid regions, groundwater resources are prone to depletion due to excessive water use and little recharge potential. Especially in sand dune areas, groundwater recharge is highly dependent on vadose zone properties and corresponding water fluxes. Nevertheless, vadose zone water fluxes under arid conditions are hard to determine owing to, among other reasons, deep vadose zones with generally low fluxes and only sporadic high infiltration events. In this study, we present an inverse model of infiltration experiments accounting for variable saturated nonisothermal water fluxes to estimate effective hydraulic and thermal parameters of dune sands. A subsequent scenario modeling links the results of the inverse model with projections of a global climate model until 2100. The scenario modeling clearly showed the high dependency of groundwater recharge on precipitation amounts and intensities, whereas temperature increases are only of minor importance for deep infiltration. However, simulated precipitation rates are still affected by high uncertainties in the response to the hydrological input data of the climate model. Thus, higher certainty in the prediction of precipitation pattern is a major future goal for climate modeling to constrain future groundwater management strategies in arid regions.

In arid regions, groundwater resources are prone to depletion due to excessive water use and little recharge potential. Especially in sand dune areas, groundwater recharge is highly dependent on vadose zone properties and corresponding water luxes. Nevertheless, vadose zone water luxes under arid conditions are hard to determine owing to, among other reasons, deep vadose zones with generally low luxes and only sporadic high iniltration events. In this study, we present an inverse model of iniltration experiments accounting for variable saturated nonisothermal water luxes to estimate effective hydraulic and thermal parameters of dune sands. A subsequent scenario modeling links the results of the inverse model with projections of a global climate model until 2100. The scenario modeling clearly showed the high dependency of groundwater recharge on precipitation amounts and intensities, whereas temperature increases are only of minor importance for deep iniltration. However, simulated precipitation rates are still affected by high uncertainties in the response to the hydrological input data of the climate model. Thus, higher certainty in the prediction of precipitation pattern is a major future goal for climate modeling to constrain future groundwater management strategies in arid regions.
In most arid countries groundwater is the primary freshwater source since surface water resources are limited and temporally extremely variable, resulting from highly temporal luctuations 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 quantiication of iniltration 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 iniltration 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 afected 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 iniltrating water have to be considered for estimating direct recharge rates in these areas (Scanlon et al., 1997). Eforts have been made to better understand and quantify vadose zone processes and water luxes by improved measurement techniques of water content (e.g., Dahan et al., 2003;Jones et al., 2005), observation of scale efects (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 diicult ield conditions, such as highly Vadose zone water luxes in arid settings are investigated regarding their sensitivity to hydraulic soil parameters and meteorological data. The study is based on the inverse modeling of highly deined soil column experiments and subsequent scenario modeling comparing different climate projections for a deined arid region.
p. 2 of 14 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 inluence of temperature gradients as well as ambient soil water contents close to residual saturation have to be considered carefully if nonisothermal water luxes 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 low are oten measured under multistep equilibrium and isothermal conditions (Wildenschild et al., 2001) and are not appropriate under nonisothermal conditions at ield 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 iniltration 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 iniltration, 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 speciically 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). he 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 speciic 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 speciic to the vadose zone. Reasons for neglecting vadose zone processes in large-scale groundwater recharge models are mainly increased computational requirements and the diiculty of inferring the soil hydraulic parameters for unsaturated zone low calculations and the subsequent parameter estimation requirements (Vrugt et al., 2004). Currently, investigations coupling long-term regional shits in climate with process-orientated, small-scale vadose zone models are still lacking (Candela et al., 2012). his 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 luxes at the more local scale of the vadose zone. hus, not only climate and regional hydrological patterns have to be considered to study recharge in arid regions properly, but also, water, vapor, and heat luxes in the context of hydraulic and thermal soil properties must be speciied.
In this study, we present an inverse model of iniltration experiments accounting for variably saturated nonisothermal water luxes 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. he 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 efect of diferent climate change projections based on greenhouse gas (GHG) emission scenarios on groundwater recharge in nonvegetated sand dune areas in an arid environment. he inverse model estimates hydraulic and thermal parameters and assesses their sensitivity for simulating heat, water, and vapor luxes under transient arid climate conditions. he scenario modeling couples global climate model projections with regional recharge estimates. his allows insight into the inluence of changing climatic conditions, as well as to derive threshold values for iniltration, evaporation, and percolation rates. hree GHG emission scenarios (A2, A1B, and B1) (Nakićenović and Swart, 2000) that give global climate projections for 2100 were investigated. hey were simulated with the coupled atmosphere-ocean model EGMAM (Huebener et al., 2007). his allows transferring the inluence of expected greenhouse gas emissions into future changes in groundwater recharge. hus, a holistic analysis of hydraulic processes in deep vadose zones in response to shits in climate in a present-day arid setting is obtained.

Investigated Arid Setting
he 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). he 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. he 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 −1 , 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 −1 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.

Data Collection
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. he 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 efect of diferent GHG emission scenarios and corresponding predicted changes in climate on downward water luxes in the deep vadose zone of the investigation area in south-eastern Saudi Arabia.

Soil Column Experiments
Large column experiments (1-m length) were conducted with repacked homogenous red dune ine 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). hese measurements were performed on repacked 100-cm 3 samples. he samples were packed to a bulk density of 1700 kg m −3 as measured in the ield. Saturated hydraulic conductivity was measured with a constant head permeability test (Wietsma et al., 2009). he 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). he 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). he Brooks and Corey model was chosen due to its validity for coarse grained material (e.g., Schaap et al., 2001). he 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 3 cm −3 . Nevertheless, residual water content can be better deined as remaining water content at high capillary pressures. hus, another itting 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 relect particularly dry conditions better. his was also expected to provide more suitable initial hydraulic parameters for the inverse modeling. hermal conductivity was measured by a transient line heat source method with a KD2 thermal analyzer (Decagon Devices) for ive steps from dry to fully saturated conditions. he 30-mm needle sensors were installed in the 100-cm 3 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 diferences in sample temperature, inluences of the measurement device, and nonequilibrium conditions caused by difusive efects (Smits et al., 2013).
For the large column experiments, the dune sand was packed with a bulk density of 1700 kg m −3 and a water content of 0.07 m 3 m −3 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. he air-illed head space at the top of the column was used to generate a low of air, which was controlled in velocity and relative humidity applying a lowmeter and forcing the air low along a speciic saturated salt solution to constrain the humidity. By measuring the relative humidity and the temperature of the inlowing and the outlowing air, evaporation and soil water uptake were calculated. he suction plate enabled column outlow 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 proiles, as well as water luxes at the column boundaries, yielding evaporation and discharge, were measured at high temporal resolution of up to 1 min. he 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 lux 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. he 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. hese irrigation patterns were chosen (i) to have an initial percolation throughout the whole column for analyzing wetting processes and (ii) for analyzing percolation and iniltration depths of lower irrigation amounts. he temperature at the top of the column was set to 40°C during the irst 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 luxes. At the bottom of the column temperature was ixed to 18°C to mimic the neutral temperature zone. he laboratory column experiments are described in detail in Pletschinger et al. (2012).

Global Climate Model Predictions for Southeastern Saudi Arabia between 1860 and 2100
To take into account the changing climate in Saudi Arabia until the end of the century, global simulations of future climate were analyzed. he 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 (Wolf et al., 1997). In the EGMAM setup ECHAM4 supports a horizontal resolution of approximately 3.75° ´ 3.75°. he 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 reinement near the equator up to 0.5°. For coupling, a lux correction is applied for heat and freshwater exchange. he orography of the grid cells is based on a digital elevation model that is averaged Table 1. Hydraulic parameters of the dune sand obtained by itting the Brooks and Corey model (1964) to small-scale retention measurements and from the inverse model of the large column experiments. Scaled dimensionless parameter sensitivity coeicients of the estimated model parameters are calculated with PEST. q r , residual water content; q s , volumetric water content at saturation; h b and l , Brooks and Corey parameters; K s , saturated hydraulic conductivity.
Hydraulic parameters of the dune sand Diferent scenarios of future greenhouse gas emissions were developed for the IPCC hird Assessment Report, the IPCC Fourth Assessment Report, and described in Nakićenović and Swart (2000). he 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 diferences in per capita income, and the rapid introduction of new and more eicient technologies with a balance across fossil and nonfossil 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. he 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 luxes, grid boxes covering the study area have to be extracted. he 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 lowlying 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.

Governing Equations of Vadose Zone Processes
Simulations were performed with Hydrus-1D (Version 4.14) (Šimůnek et al., 2009), which had previously been applied for studying water, vapor, and heat flux in arid conditions (e.g., Walvoord and Scanlon, 2004). Hydrus employs a modified Richards equation to simulate one-dimensional water, vapor, and heat low in variably saturated porous media as given in Saito et al. (2006): where x is the spatial coordinate position (m), q is the volumetric water and vapor content (m 3 m −3 ), h is the pressure head (m), T is the temperature (K), K Lh and K LT (m s −1 ) are the isothermal and thermal hydraulic conductivities of the water phase, and K vh and K vT are the isothermal and thermal hydraulic conductivities of the vapor phase due to gradients in h and T, respectively. he thermal hydraulic conductivity accounts for two efects: (i) the temperature dependence of the soil water retention curve and (ii) the temperature dependence of the surface tension of soil water (Nimmo and Miller, 1986).
he water retention curve and the unsaturated hydraulic conductivity function of Corey (1964 and1966) were used to describe the isothermal soil water retention curve: where S e is the efective saturation (-), q r is the residual volumetric water content (m 3 m −3 ), q s is the volumetric water content at saturation (m 3 m −3 ), h b (m) is a fitting parameter that can be physically related to the air-entry pressure, and l (-) is a itting parameter that can be physically related to the pore size distribution (Hillel, 2004).
The isothermal unsaturated hydraulic conductivity is defined according to: where K s is the saturated hydraulic conductivity of the water phase (m s −1 ) and l is the tortuosity (-).
For describing thermal conductivity k (W m −1 K −1 ) as function of volumetric water content the equation of Chung and Horton (1987) was used: where b1, b2, and b3 are empirical parameters (W m −1 K −1 ).
he 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). he hydraulic and thermal equations are solved in a fully coupled way.

Inverse Modeling to Estimate Vadose Zone Soil Properties
The inverse model was set up to estimate the hydraulic soil properties from column experiments (Fig. 2). he 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.
he hydraulic boundary on top of the column accounted for hourly changes of precipitation and potential evaporation. he 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 lux 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 ixed to 18°C. A suction plate was installed at the bottom of the column, which was represented in the model by deining hydraulic properties diferent from the soil material for a 0.1-m-thick bottom layer. Hydraulic properties of the suction plate were initially deined 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. hermal 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 low. A seepage face was assigned to the column bottom to mimic the suction plate. he seepage face boundary implies that zero-lux boundary conditions prevail as long as the local pressure head at the bottom of the proile is below the ixed seepage face pressure head value. As the local bottom pressure head reaches the deined 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 speciied 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 3 m −3 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 ater packing (Pletschinger 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. he parameter optimization process in PEST uses model outputs and corresponding observation measurements by minimizing the weighted sum of squared diferences. 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 luxes. 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 q r , the saturated volumetric water content q s , the Brooks and Corey parameters h b and l (Eq.
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 ixed to the value derived from the laboratory measurements of the soil retention curve and thermal conductivities, respectively.

Prediction of Present-Day and Future Groundwater Recharge
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 deined to account for a deep groundwater surface beneath the vadose zone. he bottom temperature was ixed 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 efect of long-term arid conditions on the soil water proile. For the spin-up period, monthly varying climate data obtained from the EGMAM simulations for the period 1860 until 2000 based on the dCO 2 experiment were assigned to the top of the model. he initial water content at the beginning of the spin-up period was set to 0.06 m 3 m −3 to account for an initially wetted proile. Temperature distribution was set homogenously to 18°C.
For the scenario modeling of the vadose zone water luxes from 2000 until 2100, boundary conditions at the top of the model were deined with daily variable input data. Atmospheric boundary conditions were extracted from the EGMAM simulations and captured precipitation and mean air temperature. he 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.
he calculation of the evaporation rate E 0 (mm d −1 ) combines a radiation term E rad (mm d −1 ) and an aerodynamic term E aero (mm d −1 ) (Eq. [5]) (Šimůnek et al., 2009 where l is the latent heat of vaporization (MJ kg −1 ), R n is the net radiation at surface (MJ m −2 d −1 ), r is the atmospheric density (kg m −3 ), c p is the speciic heat of moist air (kJ kg −1 °C −1 ), (e a − e d ) is the vapor pressure deicit (kPa), e a is the saturation vapor pressure at temperature T (kPa), e d is the actual vapor pressure (kPa), r a is the aerodynamic resistance (s m −1 ), D represents the slope of the vapor pressure curve (kPa °C −1 ), and g the slope of the psychrometric constant (kPa °C −1 ) (Šimůnek et al., 2009).

Vadose Zone Hydraulic and Thermal Parameters and Sensitivity Analysis
Mean residuals for water content were 0.02 m 3 m −3 (0.63 m depth), 0.025 m 3 m −3 (0.23, 0.33, 0.43, 0.53, 0.73 m depth), 0.03 m 3 m −3 (0.13 m depth), and 0.035 m 3 m −3 (0.03 m depth). he higher residuals at the top of the column resulted from the irst measurements ater 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 luctuations 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 inal optimization runs.
For the dune sands, the highest dimensionless scaled sensitivity coeicients were calculated for air-entry pressure. he air-entry pressure deines the pressure head at which the water content drops Vadose Zone Journal p. 8 of 14 below q s and hydraulic conductivity drops below K s as air enters the pores that are formerly completely illed with water. his high sensitivity could probably be explained by the fact that the airentry pressure shits the entire retention curve and the unsaturated hydraulic conductivity function. Especially in sandy soils, water retention and hydraulic conductivity signiicantly decrease ater air enters the porous structure, indicated by the steep slope of the retention curve (Fig. 4). For the ine 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). he inverse modeling conirmed the range of the airentry 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 deines the change from fast downward water drainage to capillary driven decreasing water luxes. 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 luxes that are inluenced 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 coeicient l of the Brooks-Corey soil model. he low sensitivity of both parameters probably relects the rather dry conditions that prevailed during the experiments whereas most of the studies which detected a high sensitivity of K s and l rather focus on the wet part of the retention curve. For the suction plate, dimensionless scaled sensitivity coeicients 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 coeicients 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 3 m −3 ) 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 Regarding observation data used in the calibration process, sensitivities were highest for water content observations at the time of breakthrough of the iniltration water in the corresponding depths, whereas top observations had higher sensitivities with respect to deeper observations. During gravitational low, sensitivities declined, and increased again when luxes decreased. For these two stages, bottom observations had higher sensitivities than top observations. During the drying phase, where water luxes converged to zero, sensitivities decreased again, whereas sensitivities at the top were higher presumably due to upward evaporative luxes. herefore, if monitoring campaigns are designed ater precipitation events, shallow parts of iniltration 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 proile will also provide sensitive and valuable data sets to simulate vadose zone water luxes in an arid setting.

Groundwater Recharge with Respect to Shifts in Climate
Estimated Shifts in Climate for 2100 he 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. he total diference 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 −1 (Fig. 5). his 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 −1 . he analysis of the diferent 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 signiicant trend in precipitation. herefore, changes in the GHG concentration do not strongly force precipitation over Saudi Arabia, but mainly temperature.
he strong variation in daily as well as annual precipitation and the diferences 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. here is also no common trend identiied comparable to that of temperature. his is valid for all widely used global climate models due to uncertainties in hydrological response (Christensen et al., 2007).

Impact of Greenhouse Gas Emission on Water Resources under Arid Climate
During the spin-up period, from 1860 to 2000, the water content in the 30-m-deep proile converged to steady-state conditions with a water content of 0.047 m 3 m −3 which deceased to lower water contents down to 0.02 m 3 m −3 in the upper 2 m. he water content below 2 m relects ield capacity conditions, where water is held against gravitational potential. In the upper 2 m, water content is additionally inluenced by the atmospheric conditions. Here, depending on precipitation amounts and atmospheric pressures that inluence evaporation, water luxes are either directed down-or upward. At the end of the spin-up period, the discharge reached steady-state conditions with 13 mm yr −1 . 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 proile and mean recharge rates varied around 22 to 29 mm yr −1 . Soil water storage and bottom lux, and thus recharge rates, are highly afected by precipitation pattern. he 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. he highest recharge peak was observed for the B1 scenario with 55 mm yr −1 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 −1 were also observed for the B1 scenario. For the A1B scenario, maximum and minimum recharge rates were 51 and 13 mm yr −1 in the year 2045 and 2030, respectively. In the A2 scenario, maximum recharge was obtained in 2075 with 51 mm yr −1 , 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 −1 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 diferent 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 proile for all three scenarios. The absolute value of soil water storage gradually increases from around 1.4 m 3 to around 1.55 m 3 in the whole soil proile from 2000 to 2100. his is contributed by the daily precipitation values that account for lower evaporation losses due to deeper iniltration used in the predictive modeling from 2000 to 2100, compared to averaged monthly precipitation values assigned during the spin-up period. his 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 proiles for all three scenarios basically show similar percolation patterns as a response to high precipitation events, even though the timing of the single events difers. 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 iniltration 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). Fig. 5. Simulated annual sum of precipitation and averaged temperature for the three greenhouse gas emission scenarios A1B, A2, and B1 until 2100 obtained with the climate model EGMAM. he realization used for the predictive vadose zone modeling is given in continuous lines (temperature) and solid bars (precipitation), dashed lines (temperature) and empty bars (precipitation) show two additional climate model realizations.

p. 11 of 14
When zooming in to the upper 5 m of the soil proiles, 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 signiicant increases in evaporation depth at the end of the simulation period (Fig. 8).
Comparing water and vapor luxes in summer and winter show diferent lux patterns driven by temperature (Fig. 9). In summer time, in the deeper part of the vadose zone, both water and vapor luxes are mainly directed downward. Only in the upper few centimeters of the soil proile, vapor lux is directed upward, contributing to evaporation. During winter, the vapor lux is directed upward down to 1 m soil depth, which is driven by thermal vapor luxes, as the thermal gradient is negative and thus directed upward. Water luxes during winter are directed downward in the whole soil proile, as precipitation mainly occurs in the winter time leading to gravitationally driven downward water low. Below the irst centimeter of the soil proile, both summer and winter water luxes are directed downward. In the deeper part of the soil proile (blow 2 m depth) summer water luxes approach zero, while winter water luxes reach steady-state conditions. hus, only during winter month iniltrated precipitation will reach deeper parts of the vadose zone in higher amounts and thus can enable groundwater recharge. 6 Conclusions he inverse modeling based on sensitivity coeicients revealed that highly relevant parameters for groundwater studies in arid regions are air-entry pressure and residual water content, rather than saturated hydraulic conductivity. hus, measurements of the retention curve near residual saturation can be more helpful for characterizing water low patterns in arid regions than for conducting iniltration tests. Nevertheless, measuring retention curves in the small columns still might not relect the retention curve valid under ield conditions. Transferring the results obtained by soil column experiments to ield scale still lacks in taking heterogeneities, as well as other environmental factors like vegetation, for Fig. 6. Simulated annual recharge, cumulated recharge, and changes in soil water storage with respect to the three greenhouse gas scenarios A1B, A2, and B2. Fig. 7. Simulated contour plots for water content changes and corresponding recharge over the whole simulation time for the A2 scenario realization.
p. 12 of 14 example, into account. hus, 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 −1 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 inluences 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 afect 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.
he study showed that groundwater recharge is mainly governed by precipitation rates, in which single high precipitation events, which can deeply iniltrate, are the main source for direct recharge in an arid environment. Diferent climate scenario studies corresponding to different GHG emission scenarios showed only marginal efects on recharge rates, and thus cumulative groundwater recharge rates until 2100 vary only minor for diferent climate predictions.