Quantifying the uncertainty in future groundwater recharge simulations from regional climate models

This study aims to show how future groundwater recharge (GR) simulations in arid areas respond to uncertainty in climatic parameters—a question, if explored, that bridges a gap in water resources management plans. To this aim, eight regional climate models (RCMs) under two representative concentration pathways (RCP4.5 and RCP8.5) projected four climatic parameters [surface air temperature, precipitation, wind speed, and potential evapotranspiration (PET)] over Qatar during the period of 2071–2100. Using topographic and groundwater data, a physically based water balance model was built to simulate future GR under these 16 scenarios. Results show high uncertainty in climatic parameters. Relative to the reference period (1976–2005), values varied under RCP4.5 (RCP8.5) from +1.8 to +3.4 (+3.8 to +5.6)°C for average temperature, −48% to +15% (−60% to +6%) for annual precipitation, −0.23 to +0.1 (−0.27 to +0.04) m/hour for wind speed, and from −5.7 to +12.8 (+4.3 to +17) mm for annual PET. Uncertainty in climatic parameters caused great uncertainty in future GR estimations. During the late 21st century, GR simulations varied from −67% to +64% with an average value of −20% under RCP4.5, and from −81% to +8% with an average value of −36% under RCP8.5. The greatest uncertainty resulted from the driving model, whereas the choice of emission scenario had a secondary impact. Since GR is a critical component of feeding arid aquifers, the study's findings emphasize the importance of both considering the uncertainty associated with climatic parameters and the regional climatic information chosen.


| INTRODUCTION
Minimal rainfall and lack of surface water continue to be significant concerns in arid areas, often making aquifers the primary source of potable and agricultural water. Therefore, it is vital, to manage aquifers to achieve socioeconomic and ecosystem sustainability. This management becomes a top priority if aquifers are put under immense pressure by climate change and anthropogenic activities such as population growth, rapid urbanization, and overexploitation. However, aquifers' sustainability cannot be achieved without a precise estimation of groundwater recharge (GR). Therefore, an accurate assessment of GR under varying climatic and environmental conditions is necessary for water resources management, climate change mitigation, and adaptation plans (Anibas et al., 2016;Pozdniakov et al., 2020;Xie et al., 2018). In contrast, unreliable GR assessment may make integrated water resource management impossible, leading to a significant deficit in water requirements, disruption in soil conditions, and land degradation (Moiwo et al., 2010;Xie et al., 2018).
Direct measurements (e.g., lysimeters, tensiometers, and time domain reflectometry probes) are optimum ways to calculate accurate GR in arid areas. However, these measurements are elaborate and expensive, which limits their accessibility in several countries. Additionally, these measurements are only locally applicable. Some arid areas have karst aquifers. Karst aquifers have their unique hydrogeological composition (e.g., heterogeneous permeability of aquifer matrix and conduits, heterogeneous hydraulic conductivity), making extrapolating local and wide range GR estimations untrustworthy (Martos-Rosillo et al., 2015). Further, arid environments mostly get heterogeneous, intense, short-duration rainfall (Ajjur & Al-Ghamdi, 2021a;IPCC, 2021). Such rainfall characteristics hinder determining soil moisture capacity, making the runoff and infiltration rain distributions uncertain. To overcome these issues, some other techniques are used for regional GR evaluation. These techniques include soil water balance, energy water balance, chloride mass balance, and inverse calibration in numerical modelling. All these techniques generally depend on climate, surface, and groundwater components as inputs. Uncertainty in these components can lead to significant uncertainty in GR estimations (Kurylyk & MacQuarrie, 2013;Martos-Rosillo et al., 2015).
Previous literature in arid areas focused on GR uncertainty associated with surface and groundwater components. Xie et al. (2018) randomly sampled vegetation parameters in semi-arid southeast Australia, generating 10 000 realizations of long-term GR using energy and water balance modelling. The resulted annual GRs varied between 7 and 144 mm despite all these realizations having the same climatic inputs. Martos-Rosillo et al. (2015) reviewed GR assessment in 51 carbonate aquifers in southern Spain. They found six estimation methods in semi-arid regions to have significant uncertainty. Both studies (Martos-Rosillo et al., 2015;Xie et al., 2018) provided no evidence on how GRs simulations will change in response to climatic uncertainty. Jackson et al. (2011)  They found that precipitation uncertainty in GCMs could considerably change future GR simulations. Though, recent improvements in the climate arena have advanced regional climate downscaling techniques, producing regional climate models (RCMs) that are supposed to reduce the uncertainty in climatic projections, leading to more accurate future GR simulations. Little effort has been devoted to defending this statement. Therefore, there is a pressing need to quantify the uncertainty in future GR simulations resulting from RCMs projections.
This study is comprised of two objectives: (1) quantifying the uncertainty that GR projections in arid areas are exposed to due to changes in climatic components, specifically the changes in precipitation, air temperature, wind speed, and potential evapotranspiration (PET); and (2) estimating future GRs in Qatar during the late 21 st century. To this end, the study obtained climatic projections from eight RCMs participating in the COordinated Regional climate Downscaling EXperiment (CORDEX) under two Representative Concentration Pathways (RCP4.5 and RCP8.5). Accompanied by topographic and groundwater data, these components were utilized in a physically based water balance model (WetSpass) to simulate GR during historical and future periods. The period between 1976 and 2005 was considered a reference, while 2071-2100 represented the late 21 st century.

| CLIMATE AND HYDROGEOLOGICAL DESCRIPTION
Qatar is a hyper-arid country covering 11 500 km 2 in the eastern Arabian Peninsula (Figure 1). Qatar's surface is almost flat, with an elevation that varies between 0 and 92 m above mean sea level (AMSL).

| Historical climate
Qatar has an extremely long, hot summer and a mild winter (Ajjur & Al-Ghamdi, 2021a, 2021b. The trends of average, maximum, and minimum air temperatures, precipitation, and wind speeds from 1976 to 2005 are depicted in Figure S1. The trends at six meteorological sta- there was a minor decrease in annual rainfall values (R 2 = 0.02). The wind speed trend also decreased. Figure S2 depicts the mean spatial distribution of climatic parameters over Qatar during the reference period (2071-2100). It shows evidence of significant spatial variability in climatic components, especially for maximum temperature and rainfall. For example, the mean historical maximum temperature of the six stations varied between 29.5 and 34.6 C, while the mean historical annual rainfall varied between 54.9 and 84.1 mm.

| Geological description
As Figure 1b shows, Qatar's soils are generally lithosols; they are composed of thin, calcareous, sandy loams (10-30 cm) covered by limestone debris and bedrock (Eccleston et al., 1981). Qatar's land-use patterns can be classified in three ways: bare areas, urban areas, and agricultural areas. Bare areas constitute the majority of Qatar, where there are exposed rocks and sand dunes.
Urban areas include built-up surfaces like buildings, roads, construction sites, and industrial communities. This class is mainly found in metropolitan Doha (the central-eastern part of the country) and small communities in the south and north. Finally, agricultural areas are generally spread throughout the northern part of the country.
Qatar has three principal aquifers: the northern basin, the southern basin, and the southwest basin (Eccleston et al., 1981;Schlumberger Water Services, 2009). Figure 3 depicts the location of these aquifers and groundwater levels. The primary aquifer is the northern basin since it receives higher rainfall and thus has higher recharge. The northern basin contains fresh groundwater lenses sitting on top of brackish and saline groundwater. The depth to groundwater level is up to 49 m from the ground surface, and its quality is better than that of the other aquifers (Southern and Southwest). Most farms are concentrated in the northern aquifer.
The northern aquifer includes approximately two-third of all wells, whereas the southern basin contains 28% (Schlumberger Water Services, 2009). These aquifers are karst carbonate heterogeneous (Baalousha et al., 2018;Schlumberger Water Services, 2009). During the last two decades, aquifers have been extremely overexploited due to anthropogenic and climate change stressors. The country's population has risen fivefold, and per-capita water consumption has reached an unprecedented rate (World Bank, 2020).
The main groundwater-associated problems are water table decline, storage reduction, quality degradation, and seawater intrusion. Eccleston et al. (1981) stated that GR occurred in Qatar under two circumstances: directly by intense rainfall (>10 mm/day) and indirectly by rain flow accumulating in land depressions (indirect recharge). Indirect recharge occurs five times more than intense rainfall (Schlumberger Water Services, 2009).

| Regional climate models
Future climatic simulations, including precipitation; average, minimum, and maximum air temperature; and wind speed, were obtained from eight RCMs participating in the CORDEX. The CORDEX uses regional downscaling methods to provide climate data at the regional scale. It has 14 geographic domains, and four of them include Qatar. These domains are Africa, Middle East and North Africa, South Asia, and Central Asia. This study selected all RCMs in a 0.44 resolution that comprise the r1i1p1 ensemble that runs through the year 2100 under the RCP4.5 and RCP8.5. Table 1  Further, several studies reported comparable, reliable results of the Hargreaves and Samani (1985) method to those obtained using the Penman-Monteith equation in various environments (Kurylyk & MacQuarrie, 2013;Martos-Rosillo et al., 2015). Therefore, the Hargreaves and Samani (1985) equation was used to calculate the daily PET as follows: where T, T max , and T min represent daily mean, maximum, and minimum temperature in C, respectively; R a represents the water equivalent of extraterrestrial radiation in mm/day; and the AHC is the adjusted Hargreaves coefficient, equal to 0.0023. Extraterrestrial radiation is computed from latitude and day data of the year (Allen et al., 1998). We

| Groundwater recharge modelling
This study calculates annual GR using the WetSpass model (Batelaan & De Smedt, 2007). WetSpass is a physically based, spatially distributed model that utilizes water balance at each grid to calculate GR, actual evapotranspiration (AET), and runoff. WetSpass results are validated through error maps, which constitute a part of the model outputs. Error maps measure the modelling error in seasonal (summer and winter) and annual water balance simulations. This error (±) should be zero, or very close to zero. WetSpass has been extensively tested, and its reliability has been proven in many regions worldwide (Gebreyohannes et al., 2013;Moiwo et al., 2009;Moiwo et al., 2010). WetSpass requires a set of climatic data (temperature, precipitation, wind speed, and PET), topographic data (topography, slope, soil type, and land use), and groundwater data (the depth to water table; Figure S3).
For each grid cell in October, there are four parts representing the non-uniformity of the land cover. These are vegetated, bare soil, open water, and impervious surfaces. For example, for the vegetated area, the water balance is calculated using a raster cell as follows: where P is the precipitation, S is the surface runoff, T is the transpiration, and I is the interception fraction. AET is the sum of two components (Equation 3): where E is the evaporation. The surface runoff (S) is related to rainfall characteristics (amount and intensity), the interception fraction, and the soil infiltration capacity. The interception fraction depends on the land cover and represents a constant percentage of annual rainfall.
The GR is the residual term of the water balance. In this study, WetSpass analysis was performed for two seasons.
These are summer (from April to September) and winter (from October to March). Geological and groundwater maps were prepared and resampled into a 30 m grid. This study obtained groundwater levels from 313 wells from fieldwork done by Schlumberger Water Service (see Figure 3). To incorporate the changes of storage on a seasonal basis, ground measurements of water level depths were used for each season (winter and summer). Next, mean values of seasonal F I G U R E 3 The groundwater level map for the shallow aquifers in Qatar includes the location of 313 monitored wells used in this study. Levels are presented in metres below the natural ground climatic variables during the historical period were computed by averaging monthly temperature and wind speed. Precipitation records were aggregated to get a seasonal sum then averaged to get a seasonal mean. Next, seasonal climatic maps were also resampled into a standard 30 m grid. The land use and soil type maps were connected with WetSpass model as attribute tables. The supplementary materials in the online system include these tables for the land use (Table S2), Soil type (Table S3), and runoff coefficient for different raster cells (Table S4). The reader can refer to more information about WetSpass computation, calibration, and validation in Batelaan and De Smedt (2007). After simulating seasonal GRs using the WetSpass model, the two seasons were summed to obtain yearly values. In similar, future GR simulations were conducted using inputs from 16 climatic scenarios (see Section 3.2). We created these simulations to examine how GR would change in response to differences in climatic parameters and quantify GR uncertainty. (10.8 mm) and RCM2 (11.9 mm). Other RCMs projected that future GR would not exceed 7.6 mm, with urban areas receiving negligible GR values (less than 1 mm).

| Historical groundwater recharge estimation
F I G U R E 4 The spatial distribution of annual groundwater recharge over Qatar during the historical period  To better quantify the changes in GR projections among RCMs, we computed the volume of GR using ArcGIS spatial analyst tools. Under RCP4.5, RCMs do not agree on the direction of GR change; six simulate a decline in future GR while two (RCM1 and RCM3) simulate an increase. The simulated changes in annual GR range from À24 to +23.1 Mm 3 , with an average value of À7.4 Mm 3 . In other words, the GR projections range from À67% to +64%, with an average percentage of À20%, relative to 1976-2005. The GRs simulated using RCM2 and RCM4 inputs were like the GR simulated using the RCMs ensemble.
All simulations (except RCP3) generated decreases in GR during the late 21 st century, relative to 1976-2005. Under RCP8.5, an average decline in future GR projections was estimated to be À13.2 Mm 3 (À36%). However, the differences between RCMs were larger under RCP8.5 than RCP4.5. Future GR simulations ranged from À29 Mm 3 to +2.7 Mm 3 (À81% to +8%). The GR simulated using RCM5 was the closest to that simulated by the RCMs ensemble. Interestingly, the RCM7 projected the highest decline in GR under both RCPs. F I G U R E 5 Annual GR distribution during the late 21 st century simulated using RCM inputs under RCP4.5. GR, groundwater recharge; RCM, regional climate model similar. However, this study had inconsistent findings compared with Baalousha (2016), who obtained an annual GR of 58.7 Mm 3 in 1980, and Harhash and Yousif (1985), who proposed a range between 21 and 166 Mm 3 (1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983). Although using different approaches to estimate GR may produce result variations, the variations should not be significant. Martos-Rosillo et al. (2015) reviewed GR simulations in 51 carbonate aquifers in southern Spain. They stated that semi-arid areas with <300 mm of rainfall had a GR of $38% of rainfall.

| DISCUSSION
Qatar has comparable hydrogeological characteristics: carbonate karst aquifers and receives an annual rainfall of 76 mm (see Section 2.1).
The historical simulated groundwater recharge  was 31 Mm 3 (40.7%) which is similar to Martos-Rosillo et al. (2015) findings. Therefore, the estimated historical annual GR in this study is acceptable. Regarding future estimations, it may be presumed that the future mean annual GR will decrease by À20% under RCP4.5, given that only two RCMs suggested a GR increase (see Figure 5).
Nevertheless, the trend of GR decrease is more evident under RCP8.5, making the À36% decline a more acceptable data point.
Few studies discuss the uncertainty associated with future GR simulations in arid areas, making assessing our results difficult. A comparison with research conducted in other environments shows some agreement and disagreement. Based on a single GCM, downscaled by two methods, under two emission scenarios, Holman et al. (2009) found an approximate range of GR uncertainty (between À14% and +37%) for loamy soil areas in East Anglia, UK. Jackson et al.  ]. Accordingly, it can be concluded that the effect of emission scenarios on GR simulations is less evident than the driving GCM effect.
The findings of this study agree with Crosbie et al. (2011), who attributed the most significant cause of uncertainty in GR projection in southern Australia to GCMs, whereas the choice of downscaling approach was less critical. Kurylyk and MacQuarrie (2013) reported a similar investigation in a humid climate in eastern Canada. They found that different emission scenarios produce little uncertainty in GR simulations. Kurylyk and MacQuarrie (2013) stated that the variability in future GR simulations first arises from the downscaling approach, second from the driving model, and third from the emission scenario. Holman et al. (2009) documented that more uncertainty arises from the downscaling process chosen than the emission scenario chosen.
Therefore, this study reasons that the uncertainty in GR projections, from most significant to smallest, is due to the variations in the driving GCMs and emissions scenarios.

| CONCLUSIONS
This study examined the uncertainty in future GR simulated using different climatic components. The study showed that climatic uncertainty arising from selecting different driving models and emission scenarios leads to uncertainty in GR simulations. As a result, the uncertainty in GR projections makes identifying the magnitude, and even the direction, of GR evolution challenging. The findings of this study suggest that future GR estimations are incompetent if climatic components are uncertain. Therefore, it is highly recommended that GR estimates be compared with ground truth measurements. Since, sometimes, ground truth measurements might be unavailable in arid areas (like this study), the simultaneous comparison might be impossible. Thus, there is a pressing need to enhance climate services at a regional scale, giving arid areas a vital opportunity to project reliable climatic simulations that reduce GR uncertainty and support sustainable water resources management. The focus of attention should be made on the greatest source of uncertainty, that is, driving GCMs.
Other recommendations to improve GR analysis include considering the human impact, such as changes in groundwater abstraction, soil conditions, socioeconomic aspects, and land use. Further studies should also consider seasonal and spatial variations of GR estimations.
Knowing such variations is essential to allocating abstraction projects and prioritizing augmentation schemes.
To ensure the sustainable development of Qatar aquifers, maximum abstraction should not exceed long-term recharge. Considering the projected decrease in GR, abstraction works should be minimized, and decision-makers should create and apply managed aquifer recharge projects (Ajjur & Baalousha, 2021). Additionally, the predicted decline in GR rates may imply larger runoffs under same AET losses (i.e., urban flash floods) because, in shallow aquifers, precipitation is partitioned between AET, runoff, and GR. The consequences of such results would be devastating for the built-up environment of Doha. Therefore, policymakers must implement rigid and effective regulations for effective urban planning and flood prevention.