Isotopic Characterization of Rainwater for the Development of a Local Meteoric Water Line in an Arid Climate: The Case of the Wadi Ziz Watershed (South-Eastern Morocco)

Abstract

For any hydrological or hydrogeological system, the arrival of new rains is the input signal to the system. This isotopic signature of precipitation is of major interest in understanding the recharge processes of the aquifer system. On the scale of a given basin, staged stations at different altitudes and spread out in space allow this input signal to be well characterized and to draw the local meteoric water line. In south-eastern Morocco, specifically, in the Errachidia region, several chemical and isotopic studies of the waters of the various aquifers have been carried out. In the absence of a local meteoric water line, these studies were based on the use of the global meteoric water line (GMWL). Thus, the objective of this work is the isotopic characterization and the elaboration of the local meteoric water line of the rainwater of the Ziz watershed. This characterization of the input signal in the study area is based on 41 measurements of stable isotopes (δ¹⁸O and δ²H) relating to the precipitations collected during the period from December 2019 to November 2020 in four staged stations at different altitudes and spread over the space from upstream to downstream of the watershed. The linear relationship of δ²H as a function of δ¹⁸O describes the local meteoric water line (LMWL) by equation δ²H = 7.5 ± 0.3 δ¹⁸O + 4.6 ± 1.7; R² = 0.93. This equation displays evaporation confirmed by the arrival of continental currents in an arid environment. The variation in precipitation δ¹⁸O as a function of the sampling altitudes for the rains highlighted the relationship δ¹⁸O = −0.0026 ∗ Z − 1.67, with R² = 0.93, which means an altitudinal gradient of −0.26‰ per 100 m of altitude. In this regard, the development of the local meteoric water line and the determination of the altitudinal gradient for the first time in this arid to semi-arid region of the watershed will be of great use to researchers and water resource managers; for example, to help determine the groundwater recharge areas, determine the exchanges between surface water and groundwater, and analyze many other hydrological problems.

1. Introduction

For at least four decades in the Maghreb, climate change has resulted in a gradual rise in temperatures and a decrease in precipitous waters [1,2,3,4,5,6]. In sub-arid climatic zones, establishing a precise hydrological balance capable of allowing the simulation of scenarios of active management of water resources is a delicate operation. Indeed, the irregularity of the temporal distribution of precipitation, as well as that of the measured water heights, prevent the use of averages in forecasts. In addition, the low density of meteorological stations, in particular in mountainous areas, which nevertheless contribute predominantly to effective rainfall, makes the spatialization of the balance unrepresentative. It is therefore necessary to resort to more global methods, depending less on punctual measurements. Isotope hydrology using chemical and isotopic data is a means of accessing the global balance and understanding of the hydrological behavior of aquifer systems [7]. These isotopic techniques have a wide field of application. They allow tracing the origin of the water; in particular, the determination of the average altitude of the recharge zone [8,9,10].

These techniques also enable determining mixtures between aquifers [11], or between surface water and groundwater [12]. On the other hand, they make it possible to determine the average age of water [13], with the use of tritium, suitable for water that has been recharged at least in part during the last decades [10,14]. In addition, it is undoubtedly in the management of water resources that isotopic tracing is best used, especially in arid and semi-arid zones, where other methods are most insufficient [15,16,17]. Furthermore, the isotopic signature of precipitation is of major interest for understanding the recharge processes of an aquifer system. Based on global precipitation, Craig (1961) [18] plotted the linear relationship of δ²H as a function of δ¹⁸O and defined the global meteoric water line (GMWL) of equation δ²H = 8δ¹⁸O + 10. This line is always considered as a first approximation as the norm and it is used for the tracing of groundwater, in the absence of local meteoric water lines. However, local meteoric water lines, based on the relationship δ²H as a function of δ¹⁸O of local precipitation waters, are very useful references for understanding and tracing the origin of local groundwater and its movements. They also reflect climatic variations and the effects of evaporation [7,19]. Deviations between the GMWL and LMWL are generally related to the varying conditions of the steam sources (deuterium excess) and to evaporation (slope of the straight line) [7]. On the other hand, the isotopic signal of rainwater depends strongly on the history of air masses, the meteorological conditions they experienced at the origin, during their journey and during the condensation of water vapor, re-evaporation of raindrops under the cloud base in unsaturated air, and participation of the recycled vapor. These geographic and temporal variations are the result of the mass-dependent isotope fractionation effect that occurs during phase changes of water bodies complicated by the systematics of the hydrologic cycle, including evaporation, condensation, transpiration, and the mixing of air masses [19,20,21,22,23].

In Morocco, several hydrogeological studies using isotopic techniques have been conducted [24,25,26,27,28,29,30]. These works have treated, among others, the origin and residence time of groundwater in the Tadla aquifer [27]; the flow velocity and age of water in the Tadla aquifer [28]; the recharge of aquifers by the landforms [1,29]; an isotopic study of the relationship between surface water and groundwater under a semi-arid climate: case of the Souss-Massa basin [24]; and the origin of geothermal waters in Morocco using multiple isotopic tracers [30]. These and other studies were based on the use of the GMWL established by Craig [18]. Despite Craig’s global meteoric water line being a global average of many local meteoric water lines, a characteristic LMWL of any given locality may be quite different from the global line, as each LMWL represents the site-specific long-term covariation of hydrogen and oxygen stable isotope ratios related to local climatic factors, including the origin of the water vapor, secondary evaporation during rainfall, the seasonality of precipitation, and mixing [31]. These local factors affect both the deuterium excess and the slope. Furthermore, for any detailed regional or local study of groundwater recharge using δ¹⁸O and δ²H that compare isotopic ratios in groundwater or surface water with precipitation at specific locations, it is important to define as best as possible the LMWL. In reality, LMWLs have practical utility as benchmarks for evaluating the hydroclimatic processes and for interpreting the isotope ratios measured in terrestrial (e.g., ground, river, or lake) and biologically derived (e.g., stem and leaf) waters [32]. However, in some cases, due to their rigorously monitoring precipitation over a representative period (at least year), meteoric water lines are borrowed from the closest available monitoring station. Defining the LMWL for precipitation is an important part of groundwater investigations.

In south-eastern Morocco, specifically, in the Errachidia region, several chemical and isotopic studies of the various aquifers waters have been carried out [1,2,3,4,5,6]. In the absence of a local meteoric water line in the region of the Wadi Ziz watershed, these studies were based on the use of the GMWL (Craig 1961) [18]. Hence, the objective of the present work is the isotopic characterization and the elaboration of a local meteoric water line throughout the Wadi Ziz watershed, and consequently the definition of the linear correlation of δ²H as a function of δ¹⁸O. This developed relationship, will be a continuity of LMWL defined in the High Atlas Mountains and allows to compare these local meteoric water lines drawn in this arid to semi-arid area, from the West in the Tensift watershed [33] to the Ziz basin located in the East through the Draa basin in the center by Cappy (2006) [29]. The establishment of this relationship, in this arid region (Wadi Ziz watershed), will provide a baseline for comparisons in future stable isotope studies for this region and will be of great use to researchers and managers of water resources. It will make it possible to determine the recharge areas of groundwater, to determine the exchanges between surface and groundwater, and to analyze many other hydrological problems.

2. Materials and Methods

2.1. Study Area

This study was carried out throughout the Wadi Ziz watershed, located in the south-east of Morocco (Figure 1). It is a very large watershed with an area of 14,415.8 km². It is bounded to the North by the Moulouya watershed, to the West by the Rheris watershed, and to the East by the Guir watershed. According to the 2014 census, the region of Errachidia has inhabitants. The agricultural sector is the main source of income and employs the majority of the working population, thus ranking first in the development of the local economy. The Wadi Ziz and its tributaries (Figure 1) represent the hydrographic network. The latter is rich by its important valleys but with little development of perennial watercourses. The main collector includes:

• -. An upper watercourse, East–West direction over 122 km, becoming North–South at Kerrando. This course drains the upstream Ziz watershed (Upper Ziz), which is limited to the South by the Hassan Addakhil dam with a capacity of 380 million m³. The High Ziz area limited to the South by the great South Atlas accident is marked by significant exchanges between surface and groundwater [1,3,4,5,6].

• -. An average course, generally North–South, collects the water from the intermediate watershed (Middle Ziz) from the dam to the Erfoud raft.

• -. A lower course, generally North–South, crosses the Tafilalet plain and gets lost in the Sahara (Lower Ziz).

The rainfall regime is dependent on orographic disturbances responsible for the stormy nature in summer and ocean disturbances at the start of the winter and spring rains. The scarcity of precipitation is mainly due to the Atlas reliefs, which constitute a barrier to oceanic influences. The presence of the Atlas barrier, culminating at altitudes above 3000 m, and the intrusion of hot winds of Saharan origin are at the origin of the harshness of the climate. The opening of the region to the South accentuates aridity and evaporation phenomena [34]. The dry period often lasts up to eight months, with maximum temperatures (40 °C) achieved during the months of June, July, and August. It has a relatively wet and very cold winter, with minimum temperatures around 6 °C in January.

2.2. Regional Geology

The regional geology displays three morphostructural units forming the landscape of the region from North to South (Figure 2): the Atlasian formations, the Cretaceous basin of Errachidia, and the formations of the Tafilalet plain [35]. Atlas formations are characterized by a Jurassic-type morphology [36]. The structure is folded; it comprises folds oriented SW–NE and E–W, with faulted calcareous anticlinal axes and flat-bottomed marly synclinal basins, forming parallel valleys [37]. The lithology is varied. It has two limestone complexes likely to be good aquifer reservoirs [37]. The lower Liassic is calcareous and dolomitic. This massive or bedded mass is fractured and cracked; its thickness varies from 1400 m in the North to 700 m in the South [38,39]. The Dogger, made up of marl and limestone, potentially forms a very powerful multilayer aquifer. These two rock types are separated by the permeable layer of the Toarcian. The Cretaceous basin is made up of sandstone, sandy clay, and gypsum of the Infracenomanian. The fractured and karstified limestones of the Turonian are topped by the Senonian, which consists of alternating sand, sandstone, and marl [1,2]. To the south of the Cretaceous basin, the Tafilalet plain is a depression formed by Quaternary alluvial sedimentary deposits resting directly on the Paleozoic [35].

2.3. Rain Sampling

At the scale of a given watershed, the stations staged at different altitudes and spread out in space allow this input signal to be well characterized. Thus, the rain collection for the analyses of δ¹⁸O and δ²H was carried out at four climatic stations between December 2019 and November 2020. These stations (Figure 1) are well distributed in the Ziz watershed from upstream to downstream, at different altitudes: Zaouïa Sidi Hamza (ZSH: 1740 m), Rich (RCH: 1260 m), Errachidia (ERRA: 1050 m), and Erfoud (ERF: 840 m) (Figure 1). At each station, a rainwater collection device has been installed. This device was produced in collaboration with the technicians of the National Center of Energy, Sciences and Nuclear Techniques (CNESTEN). The sampler was designed according to the recommendations of the IAEA guide [40]. It consists of a support with a funnel, and a one-liter bottle with a 7 m hose for degassing the air (Figure 3 and Figure 4). This system reduces the contact between the atmosphere and the air inside the bottle and thus avoids too high evaporation of the sample taken. Therefore, with the help of people living near the weather stations, rainwater was collected after each event when possible. The samples were stored in double-capped polyethylene bottles previously rinsed with distilled and sampled water.

2.4. Isotope Analysis

The analyses of the oxygen and hydrogen isotopic composition were carried out at the National Center for Energy, Sciences and Nuclear Techniques (CNESTEN), with a Thermo-Fisher Delta Plus mass spectrometer (IRMS), using the well-established equilibration technique. Each sample vessel was loaded with 3 mL of water and equilibrated with respect to ultrapure reference gasses (H₂ for hydrogen and CO₂ for oxygen) in an automated equilibration device kept at 18 °C [41] and coupled to the mass spectrometer. Normalization was based on conversion and reporting of raw data to that of two laboratory standard waters analyzed with unknowns that were previously calibrated using international reference materials: Vienna Standard Mean Ocean Water 2 (VSMOW2) and Standard Light Antarctic Precipitation (SLAP2). To assure the accuracy of the isotopic values, all samples were analyzed in duplicate on different days and samples with unacceptable replicate tolerances were analyzed until acceptable statistics were obtained. A third laboratory standard was used to monitor and track the long-term quality control and laboratory performance. The ¹⁸O/¹⁶O and ²H/¹H ratios were expressed in terms of departure from the international standard for water: V-SMOW. The instrumental precision is ±0.1‰ for δ¹⁸O and ±1‰ for δ²H.

2.5. Calculation of Weighted Means of Precipitation Isotopic Composition

The precipitation-weighted regression for determining meteoric water lines was introduced by Hughes and Crawford (2012) [42]. The reasoning was that smaller amounts of precipitation are more likely to have a lower d-excess. The weighted means of the isotopic composition of precipitation from the four measuring stations were calculated by the following formula.

$\overline{{\mathsf{\delta }}^{18}\mathrm{O}}=\frac{{{\displaystyle \sum }}_{i=1}^{i=n}{\mathsf{\delta }}^{18}{\mathrm{O}}_{\mathrm{i}}\times P_i}{{{\displaystyle \sum }}_{i=1}^{i=n}{P}_i}$

$\overline{{\mathsf{\delta }}^2\mathrm{H}}=\frac{{{\displaystyle \sum }}_{i=1}^{i=n}{\mathsf{\delta }}^2{\mathrm{H}}_{\mathrm{i}}\times P_i}{{{\displaystyle \sum }}_{i=1}^{i=n}{P}_i}$

where P = precipitation; and $\overline{{\mathsf{\delta }}^{18}\mathrm{O}}$ and $\overline{{\mathsf{\delta }}^2\mathrm{H}}$: weighted means of the isotopic composition of precipitation.

2.6. Global Meteoric Line and Deuterium Excess

It is now well established [43], theoretically and on the basis of global measurements, that for precipitation not subjected to evaporation, δ¹⁸O and δ²H maintain a linear correlation close to the GMWL of Craig (1961): δ²H = 8 ∗ δ¹⁸O + 10. The recurrence of a value close to 8 for slope during the entire precipitation is due to the fact that condensation occurs at equilibrium conditions (saturation) [44]. The amount d, referred to as “deuterium excess”, corresponds to the original offset, due to the fact that the vapor comes from water previously evaporated on the ocean. This difference has a value close to +10 for stations subjected to a regime of simple oceanic disturbances. It increases when a significant proportion of the condensing vapor comes from a closed sea or from continental evaporation. This is the case, for example, for the western Mediterranean, where d is of the order of 14 ‰, reflecting the mixed origin of precipitation—the Atlantic and Mediterranean. Hence, this deuterium excess is a marker for the rain’s origin. It is calculated according to the relation d = δ²H − 8 ∗ δ¹⁸O [19].

2.7. Observation of Air Masses Trajectories

In this study, in order to observe the trajectories of air masses for 96 h before their arrival at the study area, the Hybrid Single Particle Lagrangian Integrated Trajectory model (https://www.ready.noaa.gov/HYSPLIT.php (accessed on 5 June 2021)) [45,46] was used to generate the paths traversed by the air masses that give rise to rainfall in the ZIZ basin. The position of the air mass was calculated every six hours and at three different altitudes (500 m, 1500 m, and 3000 m NGF) during the 96 h preceding the rain collection. Thus, three main sectors are highlighted during the four rainy episodes observed during the study period: 6 December 2019, 22 January 2020, 19 March 2020 and 1 September 2020.

3. Results

3.1. Origin of Air Masses

The three sectors identified are as follows (Figure 5):

• -. ANW (Northwest Atlantic): this sector concerns the air masses coming from the Northwest Atlantic Ocean and crossing the land surface, the Southeast of Morocco, before reaching the sampling area.

• -. EAST Europe: These air masses move over Eastern Europe, the Mediterranean Sea, Tunisia, and Algeria.

• -. S-SE (south and south-east): These air masses come from the south-east of the area. They therefore pass over Algeria and pass over south-eastern Morocco before reaching the sampling area.

The analysis of the retro-trajectories of the air masses show that the majority of precipitation in the ZIZ watershed comes from the Atlantic, Eastern Europe, and the Mediterranean Sea, and are mainly influenced by continental circulations in the east and south-east of Morocco.

3.2. Climatic Conditions at the Sampling Stations

The rains in the Ziz watershed are generally short and brutal. They last no more than two to three days without stopping. The total amount of rain measured over different time durations in the study stations is 251.7 mm at the Zaouïa Sidi Hamza station with a number of 45 rainy days against 113.4 mm and 26 rainy days in Erfoud (Table 1 and Figure 6 and Figure 7)). Between these two stations, the Errachidia station located in the intermediate watershed received a quantity of rain of 122.4 mm with 31 rainy days and the Rich station upstream in the High Atlas with 202.4 mm and 58 days of rain (Table 1 and Figure 6 and Figure 7).

The isotopic composition of precipitation depends on the different climatic conditions that prevail in the different sub-watersheds of the Wadi Ziz. Indeed, a temperature difference and humidity gradients exist between the highest station located to the North in the High Ziz sub-watershed (ZSH: 1760 m, mean T = 15.95 °C mean air moisture = 48%) and the Erfoud station located in the Lower Ziz to the south (840 m, T = 21.98 °C and mean air moisture = 24%). High variability in climatic conditions along the altitudinal gradient have also an impact on the isotopic composition of precipitation at different elevations. Between these two stations, the Rich station (1240 m, average yearly T = 17.98 °C and mean air moisture = 24%) located in the upper Ziz and the Errachidia station (1050 m, average yearly T = 20.94 °C and mean air moisture = 36%) located in the Middle Ziz sub-watershed show intermediate values (Table 1 and Table 2; Figure 8 and Figure 9). Based on the available data range from 840 m to 1800 m of the four points in Figure 8, the mathematical expression of the relationship between annual rainfall and elevation is summarized by a straight-line with the equation of y = 0.16* × (−22.85), R² = 0.90.

3.3. Isotopic Characterization of Rainwater in the Study Area

The results of the analyses and calculation of d (deuterium excess) are shown in Table 3. The analysis of the data displays that the stable isotope contents vary between −10.27 ‰ and −1.98‰ for δ¹⁸O and between −77.50‰ and −11.70‰ for δ²H, with means at −5.09‰ and −33.70‰, respectively. The deuterium excess in water calculated according to the relationship (d = δ²H–8 ∗ δ¹⁸O) varies from 1.48 to 15.72 ‰, respectively. The mean of “d” is 7.05‰ (Table 4). Taking into account that for the northern hemisphere, the d-excess values for precipitation vary from about 5 to 15‰ [47], our values appear quite normal, while low values <5‰ indicate partial evaporation of water either during precipitation from the cloud or during sample collection and storage during warm and dry conditions.

4. Discussion

4.1. Development of the Local Meteoric Water Line (LMWL)

The graphical representation of the isotopic contents of all the available precipitations in the four stations of the study area, shown in the Figure 10, shows a very good distribution of the measurements data and exhibit a tight linear relationship between δ²H and δ¹⁸O. Hence, with the objective of determining a consistent LMWL to be used for a variety of hydrological applications in this continental region, a precipitation weighted regression method [42,48] aiming to reduce the effect of rainfall amount of small precipitation samples on the calculation of the LMWL was used for the development of an updated LMWL with the equation δ²H = 7.5 ± 0.3 ∗ δ¹⁸O + 4.6 ± 1.7 and a correlation coefficient of R² = 0.93. The obtained relationship using the ordinary least squares regression (OLSR), is similar to that established in the Drâa basin by Cappy in 2006 [29], δ²H = 7.12 ∗ δ¹⁸O + 1.54, using the same method.

The slope of 7.5 ± 0.3 and intercept of 4.6 ± 1.7 of the LMWL are lower than those of the GMWL, indicating that precipitation are affected by non-equilibrium processes during formation of the atmospheric vapor masses and precipitation generating a variable difference between δ¹⁸O and δ²H because of the different kinetic effect during the fractionation of these two isotopes. Other similar situations to the produced LMWL are found in other areas of Morocco [26,28,29,49,50], revealing maybe the influence of the Atlantic, Mediterranean, and the continental disturbances characterized by a different climate, precipitation regime, and moisture sources on the whole country of Morocco.

Comparison of the data measurements with the GMWL shows that the most depleted rainwater samples are plotted very close to this line and the majority of the most enriched samples are located below it, reflecting that the sub cloud evaporation of lighter molecules from rain drops in the arid climatic conditions of the study area was an important factor affecting the stable isotopic composition of liquid precipitation falling through dry air [51]. The projection of the weighted mean values of the isotopic composition of the precipitations of the various stations indicates clearly that the rains collected at the stations of Erfoud and Errachidia are the most enriched in heavy isotopes (ERF: δ¹⁸O = −4.06‰ and δ²H = −24.26‰; ERRA: δ¹⁸O = −4.07‰ and δ²H = −26.97‰), compared to the rainfall collected at the stations of Rich and Zaouia Sidi Hamza. This is due to the low altitudes of this area of south-eastern Morocco, characterized by high aridity, which recorded the lowest humidity of the stations studied, with values of 24% and 38%, respectively. Furthermore, the arrival of continental circulations from the east and south-east of Morocco may also cause significant evaporation and thereby result in enrichment in heavy isotopes.

4.2. Seasonal Variation

Seasonal variations in isotopic signatures of precipitation are described for the four stations during the two climatological periods (Figure 11): wet (November–April) and dry (May–October). The ZSH and Rich stations experience a notable enrichment in heavy isotopes compared to the wet period (Table 5); this results in summer precipitation with more of the recycled water condensed at a warmer temperature.

4.3. Local Altimetric Gradient of δ¹⁸O as a Function of Altitude

To assess groundwater resources and to define the boundaries of aquifer protection zones, it is essential to locate the recharge zone. The altitude effect of stable isotopes often solves this type of problem in regions with marked topography. Thus, it is possible to determine the average altitude of the water supply area of a source by locating its average isotopic composition on a regional correlation line between the isotope composition of the oxygen or hydrogen and the altitude. This straight line can be defined by the relation of the ¹⁸O of precipitation as a function of the altitude of the rainfall. In fact, the isotopic composition of precipitation varies, depending on the altitude of the sampling. Precipitation becomes increasingly depleted of ¹⁸O and ²H as altitude increases [31]. Throughout the Ziz watershed, the variation in ¹⁸O of precipitation as a function of the altitudes of sampling for the rains revealed the following relationship (Figure 12): δ¹⁸O = −0.0026 ∗ Z − 1.67 with R² = 0.93, which means an altitudinal gradient of −0.26‰ per 100 m of altitude. This same value (−0.26‰) of the gradient was found by Abourida (2007) [52] in the Haouz basin and the value of −0.27‰ per 100 m was established by El Ouali et al. (1999) [1] in the Errachidia Basin. It is also close to that established by Cappy (2006) [29] in the Daraa watershed (−0.29‰ in δ¹⁸O per 100 m) for altitude and that established by Marcé (1975) [41] in the Atlas Range in Morocco, which is −0.3‰ in δ¹⁸O per 100 m. Finally, this value of −0.26‰ is slightly lower that of −0.32‰ in δ¹⁸O per 100 m established by Raibi et al. (2006) [33] in the Tensift watershed and −0.35‰ per 100 m by El Ghali et al. (2018) [50] in the Berrechid basin.

5. Conclusions

The analysis of the retro-trajectories of the air masses shows that the majority of precipitated rains in the Ziz watershed come mainly from the Atlantic, Eastern Europe, and the Mediterranean Sea, and are mainly influenced by continental circulations in the east and south-east of Morocco. The δ²H versus δ¹⁸O relationship of the rain collected from four stations characterized by heterogeneous elevation, and thus are well distributed along an altitudinal gradient, made it possible to define the input signal of the groundwater tables in the Ziz watershed, and to establish the LMWL. This line provides on the one hand a basis of comparison for future studies of stable isotopes for this region, and on the other hand, represents a major interest in the knowledge of the characteristics of recharge water over the entire watershed. Finally, an altimetric isotopic gradient estimated at −0.26‰ per 100 m of ascent in altitude is defined. Applying this gradient to the main sources makes it possible to go back to their average recharge altitudes and therefore to their hydrogeological watersheds. It is a very useful tool for the identification of recharge areas. It also allows the assessment of the hydrological processes and interactions that occur during recharge. As a perspective, this work, based on a one-year observation period, requires a longer observation period to reinforce and complete these results. Hydrochemical analysis of the rainwater is also needed.

Author Contributions

Conceptualization, A.E.O., A.R., J.M. and T.E.G.; data curation, A.E.O., A.R., T.E.G. and M.E.H.; formal analysis, A.E.O., A.R., J.M. and T.E.G.; funding acquisition, A.E.O., A.E. and A.V.R.; investigation, A.E.O.; methodology, A.E.O., A.R. and T.E.G.; project administration, A.L., M.Q., A.E.H., A.E. and A.V.R.; resources, A.E.O.; software, A.E.O. and A.R.; supervision, A.L., J.M. and A.E.H.; validation, A.E.O., A.L. and M.E.H.; visualization, A.E.O. and A.L.; writing—original draft, A.E.O.; writing—review and editing, A.E.O., J.M., T.E.G., M.E.H., A.A. and A.V.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not Applicable.

Informed Consent Statement

Not Applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors thank all those who have helped to carry out this work, which is the first in the region. We would like to thank the Ziz, Rheris and Guir hydraulic basin agency, which provided us with the climatological data; CNESTEN for the isotope analyses; and, finally, a warm thanks to the station keepers who allowed us to harvest the rains.

Conflicts of Interest

The authors declare that they have no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Geographical location and hydrographic network.

Figure 2. Geology of the study area.

Figure 3. Rain collection system.

Figure 4. Climatic station.

Figure 5. Typical retro-trajectories of air masses coming from the three defined Sectors (6 December 2019, 22 January 2020, 19 March 2020 and 1 September 2020).

Figure 6. Variation in monthly rainfall at the four sampling stations (period December 2019–November 2020).

Figure 7. Number of rainy days at the four sampling stations (period December 2019–November 2020).

Figure 8. Correlation of rain vs. altitude.

Figure 9. Annual mean values of temperature and humidity as a function of the altitudes of the four sampling stations ((ZSH = Zaouia de Sidi Hamza, RCH = Rich, ERRA = Errachidia, and ERF = Erfoud) (for their localities, see Figure 1).

Figure 10. Relationship between δ18O and δ2H of the samples collected in the four studied stations with their weighted mean values during the period December 2019–November 2020. (Zaouia Sidi Hamza = 13 values; RCH = 14 values, Errachidia = 7 values; and Erfoud = 7 values). Local Meteoric Water Line (LMWL) and the Global Meteoric Water Line (GMWL, according to Craig 1961); R2 = 0.93.

Figure 11. Seasonal variation in stable isotopes in precipitation observed during rainy episodes in the wet season (October to April) and dry period (May to September) collected between December 2019 and December 2020 at weather stations.

Figure 12. Altimetric gradient of the weighted means in δ18O according to the altitude.

Table 1

Monthly mean temperature (T °C), monthly rain (P) and number of rainy days (NRD) at the sampling stations (period December 2019–November 2020).

	Month	D	J	F	M	A	M	J	J	A	S	O	N	Total
ZSH	T °C	6.4	5.5	7.1	12.3	15.6	20.2	23.4	27.1	26.4	23.4	19.2	11.2
	P (mm)	8.8	27.6	0	36.7	10.9	25	12	26.5	13	44.6	32.2	14.4	251.7
	NRD	2	4	0	7	4	7	1	4	5	6	2	3	45
RIC	T °C	7.7	6.4	7.5	15.5	19.6	23.2	26.1	29.9	29	24.7	20.1	13.8
	P (mm)	7.4	6.5	0	49.5	12.1	17.1	3	4.9	30.2	47.9	6.4	17.4	202.4
	NRD	4	3	0	12	5	5	4	2	8	9	2	4	58
ERRA	T °C	12.6	9.39	15	17.17	20.7	25.4	29.9	34.2	34	29.1	21.2	15.2
	P (mm)	5.7	0.6	0	61.6	8.4	1.4	0.2	2.1	1	0.2	30	11.2	122.4
	NRD	5	1	0	6	3	2	1	1	3	1	3	5	31
ERF	T °C	12.8	9.1	15	17.2	20.7	25.4	29.6	33.5	33.1	29.4	22.1	15.8
	P (mm)	3.9	0.4	0	14.8	10.2	0.4	0	0	2.6	44.6	32.2	4.3	113.4
	NRD	4	1	0	4	1	1	0	0	4	4	2	5	26

Table 2

Mean climatic parameters for the period December 2019–November 2020 of the four climatic stations. (Zaouïa Sidi Hamza (ZSH), Rich (RCH), Errachidia (ERRA), and Erfoud (ERF)).

Station	Altitude (m)	Yearly Rain Value (mm)	Yearly Mean Temperature °C	Yearly Mean Relative Moisture %
ZSH	1760	251.7	15.95	48
RCH	1240	202.4	17.98	41
ERRA	1050	122.4	20.94	36
ERF	840	113.4	21.98	24

Table 3

Results of the isotope analysis, deuterium excess “d”, and precipitation (P).

	Zaouia Sidi Hamza (ZSH)				Rich (RCH)				Errachidia (ERRA)				Erfoud (ERF)
Date	δ18O	δ2H	P.mm	d	δ18O	δ2H	P.mm	d	δ18O	δ2H	P.mm	d	δ18O	δ2H	P.mm	d
2019-12-03					−7.33	−52.6	3.2	6.04
2019-12-06	−7.19	−44.2	7.8	13.32	−4.68	−21.9	3.3	15.54	−3.28	−17.3	3.9	8.94	−1.98	−11.7	2.8	4.14
2020-01-21					−3.96	−30.2	5.7	1.48
2020-01-22	−10.21	−68.6	18.7	13.08
2020-03-14	−10.02	−77.5	5.7	2.66
2020-03-18	−6.77	−50.9	3.1	3.26
2020-03-19	−6.69	−39.2	18	14.32	−3.66	−24.6	27	4.68	−5.25	−35.3	47	6.7
2020-03-23	−5.3	−29	2	13.4
2020-03-24													−5.48	−34.2	9.5	9.64
2020-03-28	−4.89	−28.7	10	10.42	−4.12	−21.8	3.6	11.16
2020-04-18					−5.26	−35.4	8.3	6.68	−5.77	−36.5	6.8	9.66	−2.98	−15.8	10.2	8.04
2020-04-24	−3.29	−24	7.5	2.32
2020-05-15	−10.27	−77.5	7.7	4.66	−5.98	−42.1	10.1	5.74
2020-05-28					−4.69	−21.8	3.5	15.72
2020-06-10	−4.35	−31.9	12	2.9
2020-07-14					−4.66	−30.1	2.6	7.18
2020-07-20									−6.02	−41.1	2.1	7.06
2020-07-21					−4.26	−29.9	2.3	4.18
2020-09-01	−4.69	−32.5	13	5.02	−5.02	−30.1	15	10.06					−5.31	−36.2	17.5	6.33
2020-09-17					−4.86	−30.1	17	8.78					−4.89	−36.8	20.5	2.32
2020-10-21									−3.96	−30.2	16	1.48
2020-10-22	−5.68	−37.7	32	7.74	−4.08	−30.6	5.1	2.04	−2.99	−17.2	12	6.72	−2.14	−14.6	24	2.52
2020-11-13					−2.14	−14.4	5.6	2.72
2020-11-19	−5.47	−38	8	5.76
2020-11-26									−4.02	−25.4	5.6	6.76	−5.28	−34.3	2.7	7.96

Table 4

Statistical parameters of the isotopic composition of the rainfall in the Ziz watershed.

Parameter	Minimum	Maximum	Mean	Standard Deviation %
δ18O, in ‰	−10.27	−1.98	−5.09	1.92
δ2H, in ‰	−77.50	−11.70	−33.70	14.94
d in ‰	1.48	15.72	7.05	3.99

Table 5

Mean values weighted with respect to the amount of precipitation of the wet and dry periods.

	Wet Period			Dry Period			Difference (Wet-Dry)
	δ18O	δ2H	d-excess	δ18O	δ2H	d-excess	δ18O	δ2H	d-excess
ZSH	−7.10	−46.56	8.73	−5.71	−39.80	5.08	1.39	6.75	−3.65
Rich	−5.08	−33.46	6.90	−3.73	−25.95	7.67	1.35	7.51	0.77
ERRA	−4.07	−26.97	8.02	−4.99	−31.79	5.09	−0.92	−4.82	−2.93
Erfoud	−4.06	−24.26	7.45	−3.94	−28.02	3.72	0.11	−3.76	−3.72

References

1. El Ouali, A.; Mudry, J.; Mania, J.; Chauve, P.; Elyamine, N.; Marzouk, M. Present Recharge of an Aquifer in a Semi-Arid Region: An Example from the Turonian Limestones of the Errachidia Basin, Morocco. Environ. Geol.; 1999; 38, pp. 171-176. [DOI: https://dx.doi.org/10.1007/s002540050413]

2. Ammary, B. Etude Géochimique ET Isotopique des Principaux Aquifères du Bassin CRétacé d’Errachidia ET de la Plaine de Tafilalet. Ph.D. Thesis; Universit Mohammed V-Agdal: Rabat, Morocco, 2007; 134.

3. Lgourna, Z.; Warner, N.; Bouchaou, L.; Boutaleb, S.; Tagma, T.; Hssaisoune, M.; Ettayfi, N.; Vengosh, A. Nitrate Contamination of Alluvial Groundwater in the Ziz Basin, Southeastern Morocco. Moroc. J. Chem.; 2014; 2, pp. 2-5.

4. Lgourna, Z.; Warner, N.; Bouchaou, L.; Boutaleb, S.; Hssaisoune, M.; Tagma, T.; Ettayfi, N.; Vengosh, A. Elucidating the Sources and Mechanisms of Groundwater Salinization in the Ziz Basin of Southeastern Morocco. Environ. Earth Sci.; 2015; 73, pp. 77-93. [DOI: https://dx.doi.org/10.1007/s12665-014-3396-1]

5. Nouayti, N.; Khattach, D.; Hilali, M. Assessment ofphysico-chemical quality of groundwaterof the Jurassic aquifers inhigh basin of Ziz (Central High Atlas, Morocco). J. Mater. Environ. Sci.; 2015; 6, pp. 1068-1081.

6. El Ouali, A.; El Hafyani, M.; Roubil, A.; Lahrach, A.; Essahlaoui, A.; Hamid, F.E.; Muzirafuti, A.; Paraforos, D.S.; Lanza, S.; Randazzo, G. Modeling and Spatiotemporal Mapping of Water Quality through Remote Sensing Techniques: A Case Study of the Hassan Addakhil Dam. Appl. Sci.; 2021; 11, 9297. [DOI: https://dx.doi.org/10.3390/app11199297]

7. Etcheverry, D. Valorisation des Méthodes Isotopiques Pour les Questions Pratiques LIéES Aux Eaux Souterraines Isotopes de l’OxygèNE ET de l’Hydrogène; FOWG Reports; Geological Series; Federal Office for Water and Geology: Berne, Switzerland, 2002; Volume 2, pp. 1-70.

8. N’da, B.; Bouchaou, L.; Reichert, B.; Hanich, L.; Danni, S.; Ait Brahim, Y.; Chehbouni, A.; Beraaouz, E.; Michelot, J.L. Chemical Tracers and Stable Isotopes Mixing Models for Groundwater Quality and Recharge Study in the Moroccan High Atlas Mountains. Groundwater and Global Change in the Western Mediterranean Area; Environmental Earth Sciences Calvache, M.L.; Duque, C.; Pulido-Velazquez, D. Springer International Publishing: Cham, Switzerland, 2018; pp. 235-244. ISBN 978-3-319-69355-2

9. Heiß, L.; Bouchaou, L.; Tadoumant, S.; Reichert, B. Multi-Tracer Approach for Assessing Complex Aquifer Systems under Arid Climate: Case Study of the River Tata Catchment in the Moroccan Anti-Atlas Mountains. Appl. Geochem.; 2020; 120, 104671. [DOI: https://dx.doi.org/10.1016/j.apgeochem.2020.104671]

10. Blavoux, B.; Letolle, R. Apports des techniques isotopiques à la connaissance des eaux souterraines. Géochronique; 1995; 54, pp. 12-15.

11. Gonfiantini, R.; Conrad, G.; Fontes, J.; Sauzay, G.; Payne, B. Etude Isotopique de La Nappe Du Continental Intercalaire et des Relations Avec les Autres Nappes Du Sahara Septentrional. Isot. Tech. Groundw. Hydrol.; 1974; 1, pp. 227-241.

12. McCarthy, K.A.; McFarland, W.D.; Wilkinson, J.M.; White, L.D. The Dynamic Relationship between Ground Water and the Columbia River: Using Deuterium and Oxygen-18 as Tracers. J. Hydrol.; 1992; 135, pp. 1-12. [DOI: https://dx.doi.org/10.1016/0022-1694(92)90078-A]

13. Fritz, P.; Fontes, J.C. Hand Book of Environmental Isotope Geochemistry; Springer: Berlin/Heidelberg, Germany, 1980; Volume 1–3.

14. Payne, B.R. The Status of Isotope Hydrology Today. J. Hydrol.; 1988; 100, pp. 207-237. [DOI: https://dx.doi.org/10.1016/0022-1694(88)90186-2]

15. Fontes, J.C.; Jean-Claude, O. Gradient isotopique entre 0 et 4000 m dans les précipitations du Mont Cameroun (Résumé de communication). Proceedings of the Quatrième Réunion Annuelle des Sciences de la Terre; Paris, France, 13–16 April 1976; Société Géologique de France: Paris, France, 1976; 171.

16. Fontes, J.C. Stable Isotope Hydrology: Deuterium and Oxygen-18 in the Water Cycle; Technical Reports Series No. 210; International Atomic Energy Agency: Vienna, Austria, 1981; pp. 273-302.

17. Fontes, J.C.; Edmunds, W.M. The Use of Environmental Isotope Techniques in Arid Zone Hydrology: A Critical Review; UNESCO: Paris, France, 1989; 75.Available online: https://snia.mop.gob.cl/repositoriodga/bitstream/handle/20.500.13000/1654/ISO729.pdf?sequence=1 (accessed on 15 January 2020).

18. Craig, H. Isotopic Variations in Meteoric Waters. Science; 1961; 133, pp. 1702-1703. [DOI: https://dx.doi.org/10.1126/science.133.3465.1702] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17814749]

19. Dansgaard, W. Stable Isotopes in Precipitation. Tellus; 1964; 16, pp. 436-468. [DOI: https://dx.doi.org/10.3402/tellusa.v16i4.8993]

20. Merlivat, L. Molecular Diffusivities of H2 160, HD160, and H2180 in Gases. J. Chem. Phys.; 1978; 69, pp. 2864-2871. [DOI: https://dx.doi.org/10.1063/1.436884]

21. Mook, W.; Rozanski, K. Environmental Isotopes in the Hydrological Cycle; IAEA Publish: Vienna, Austria, 2000; 39.

22. Araguás-Araguás, L.; Froehlich, K.; Rozanski, K. Deuterium and oxygen-18 isotope composition of precipitation and atmospheric moisture. Hydrol. Processes; 2000; 14, pp. 1341-1355. [DOI: https://dx.doi.org/10.1002/1099-1085(20000615)14:8<1341::AID-HYP983>3.0.CO;2-Z]

23. Kendall, C.; McDonnell, J. Isotope Tracers in Catchment Hydrology; Elsevier: Amsterdam, The Netherlands, 1998.

24. Bouchaou, L.; Michelot, J.L.; Vengosh, A.; Hsissou, Y.; Qurtobi, M.; Gaye, C.B.; Bullen, T.D.; Zuppi, G.M. Application of Multiple Isotopic and Geochemical Tracers for Investigation of Recharge, Salinization, and Residence Time of Water in the Souss–Massa Aquifer, Southwest of Morocco. J. Hydrol.; 2008; 352, pp. 267-287. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2008.01.022]

25. Louvat, D.; Bichara, S. Étude de Plusieurs Systèmes Aquifères du Maroc à l’Aide des Isotopes du Milieu; Rapport IAEA; IAEA Publish: Vienna, Austria, 1990; 30p.

26. Michelot, J.L. Hydrogéologie Isotopique des Systèmes Aquifères de Fès-Meknès; Errachidia et Kheng-el-Hammam–Projet MOR/8/004; IAEA Publish: Vienna, Austria, 1991.

27. Bouchaou, L.; Michelot, J.L.; Qurtobi, M.; Zine, N.; Gaye, C.B.; Aggarwal, P.K.; Marah, H.; Zerouali, A.; Taleb, H.; Vengosh, A. Origin and Residence Time of Groundwater in the Tadla Basin (Morocco) Using Multiple Isotopic and Geochemical Tools. J. Hydrol.; 2009; 379, pp. 323-338. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2009.10.019]

28. Marah, H.; Zine, N.; Qurtobi, M.; Zerouali, A. Sens d’écoulement, Vitesse et Âge des Eaux de l’aquifère Turonien du Bassin de Tadla (Maroc). Afr. J. Env. Assess Manag.; 2007; 12, pp. 1-12.

29. Cappy, S. Hydrogeological Characterization of the Upper Drâa Catchment: Morocco. Ph.D. Thesis; Universitäts und Landesbibliothek: Bonne, France, 2006.

30. Bouchaou, L.; Warner, N.R.; Tagma, T.; Hssaisoune, M.; Vengosh, A. The Origin of Geothermal Waters in Morocco: Multiple Isotope Tracers for Delineating Sources of Water-Rock Interactions. Appl. Geochem.; 2017; 84, pp. 244-253. [DOI: https://dx.doi.org/10.1016/j.apgeochem.2017.07.004]

31. Clark, I.; Fritz, P. Environmental Isotopes in Hydrogeology; CRC Press: New York, NY, USA, 1997.

32. Putman, A.L.; Fiorella, R.P.; Bowen, G.J.; Cai, Z. A Global Perspective on Local Meteoric Water Lines: Meta-analytic Insight Into Fundamental Controls and Practical Constraints. Water Resour. Res.; 2019; 55, pp. 6896-6910. [DOI: https://dx.doi.org/10.1029/2019WR025181]

33. Raibi, F.; Benkaddour, A.; Hanich, L.; Chehbouni, A.; Chtioui, M. Variation in stable isotopic contents of rainfall under semi-arid climate. The Tensift basin case study, Morocco). Proceedings of the GIRE3D—International Conference on Integrated Water Resources Management and Challenges of Sustainable Development; Marrakech, Morocco, 23–25 May 2006.

34. Kabiri, L.; Boudad, L.; Delaune, M.; Krimou, A. Nouvelles Données Chronostratigraphiques des Carbonates Continentaux Quaternaires de Tafilalet et Environs. Note et Mémoires de Géologie Marocaine; 2003; 42, pp. 341-344.

35. Margat, J. Etude Hydrogéologique du Basin Quaternaire de Tafilalet. Ressources en Eau du Maroc; 1977; pp. 310-380.

36. Joly, F. Principes pour une méthode de cartographie géomorphologique. Bulletin de l’Association de Géographes Français; 1962; 39, pp. 270-278. [DOI: https://dx.doi.org/10.3406/bagf.1962.5616]

37. Dresnay, R.D. Sédiments jurassiques du domaine des chaînes atlasiques du Maroc. Symposium sur la Sédimentation Jurassique W. Européen; Publication Spéciale Association des Sédimentologistes Français: Paris, France, 1977; Volume 1, pp. 345-365.

38. Bernasconi, R. Géologie Du Haut Atlas de Rich (Maroc). Ph.D. Thesis; Université de Neuchâtel: Neuchâtel, Switzerland, 1983.

39. Brechbuhler, Y.A. Etude Structural et Géologique Du Haut-Atlas Calcaire Entre le Jbel Ayachi et Rich (Maroc). Ph.D. Thesis; Université de Neuchâtel: Neuchâtel, Switzerland, 1984.

40. International Atomic Energy Agency (IAEA). IAEA/GNIP Precipitation Sampling Guide. 2014; Available online: http://www-naweb.iaea.org/napc/ih/documents/other/gnip_manual_v2.02_en_hq.pdf (accessed on 5 November 2019).

41. Marcé, A. Contribution des Méthodes Isotopiques à l’Étude des Modalités d’alimentation et de Renouvellement des Réserves Souterraines du Maroc; Rap. 75SGN447 French Geological Survey (BRGM): Orléans, France, 1975; Volume 131.

42. Hughes, C.E.; Crawford, J. A New Precipitation Weighted Method for Determining the Meteoric Water Line for Hydrological Applications Demonstrated Using Australian and Global GNIP Data. J. Hydrol.; 2012; 464–465, pp. 344-351. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2012.07.029]

43. UNESCOInternational Atomic Energy Agency (IAEA). Environmental Isotopes in the Hydrological Cycle: Principles and Applications. 2001; Available online: http://www-naweb.iaea.org/napc/ih/documents/global_cycle/Environmental%20Isotopes%20in%20the%20Hydrological%20Cycle%20Vol%201.pdf (accessed on 15 January 2020).

44. Gat, J.; Gonfiantini, R. Stable Isotope Hydrology: Deuterium and Oxygen-18 in the Water Cycle; International Atomic Energy Agency: Vienna, Austria, 1981; 335.

45. Draxler, R.R.; Hess, G. An Overview of the HYSPLIT_4 Modelling System for Trajectories. Aust. Meteorol. Mag.; 1998; 47, pp. 295-308.

46. Ait Brahim, Y.; Bouchaou, L.; Sifeddine, A.; Khodri, M.; Reichert, B.; Cruz, F.W. Elucidating the Climate and Topographic Controls on Stable Isotope Composition of Meteoric Waters in Morocco, Using Station-Based and Spatially-Interpolated Data. J. Hydrol.; 2016; 543, pp. 305-315. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2016.10.001]

47. Froehlich, K.; Gibson, J.J.; Aggarwal, P.K. Deuterium Excess in Precipitation and Its Climatological Significance; No. IAEA-CSP 2002-13/P International Atomic Energy Agency: Vienna, Austria, 2002.

48. Crawford, J.; Hughes, C.E.; Lykoudis, S. Alternative Least Squares Methods for Determining the Meteoric Water Line, Demonstrated Using GNIP Data. J. Hydrol.; 2014; 519, pp. 2331-2340. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2014.10.033]

49. International Atomic Energy Agency. Isotopic Composition of Precipitation in the Mediterranean Basin in Relation to Air Circulation Patterns and Climate: Final Report of a Coordinated Research Project; International Atomic Energy Agency: Vienna, Austria, 2005; ISBN 978-92-0-105305-3 Available online: https://www-pub.iaea.org/MTCD/Publications/PDF/te_1453_web.pdf (accessed on 15 January 2020).

50. El Ghali, T.; Qurtobi, M.; Marah, H.; Raibi, F.; El Mansouri, B. Chemical and isotopic characterization of rainwater in the Berrechid basin, Morocco. J. Water Sci. Environ. Technol.; 2018; 3, pp. 335-341.

51. Anuard, P.-G.; Julián, G.-T.; Hugo, J.-F.; Carlos, B.-C.; Arturo, H.-A.; Edith, O.-T.; Claudia, Á.-S. Integration of Isotopic (²H and ¹⁸O) and Geophysical Applications to Define a Groundwater Conceptual Model in Semiarid Regions. Water; 2019; 11, 488. [DOI: https://dx.doi.org/10.3390/w11030488]

52. Abourida, A. Approche Hydrogeologique de la Nappe du Haouz (Maroc) par Teledetection, Isotopie, Sig ET Modelisation; Université Cadi Ayyad Faculté des Sciences Semlalia: Marrakech, Morocco, 2007; 146p.