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. 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

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. 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

; 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. 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

. 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. 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

). 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. 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

); 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.