Next Article in Journal
Baseflow Trends for Midsize Carpathian Catchments in Poland and Slovakia in 1970–2019
Previous Article in Journal
Risk Assessment of Dike Based on Risk Chain Model and Fuzzy Influence Diagram
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Spatiotemporal Analysis of Groundwater Resources in the Saïss Aquifer, Morocco

1
Faculty of Science and Technology of Fez, Sidi Mohamed Ben Abdallah University, Fez P.O. Box 2202, Morocco
2
School of Architecture, Planning and Design, Mohammed VI Polytechnic University, Benguerrir P.O. Box 43150, Morocco
3
National School of Applied Sciences, Sidi Mohamed Ben Abdallah University, Fez P.O. Box 72, Morocco
4
Hassania School for Public Works Engineering, Casablanca P.O. Box 8108, Morocco
*
Author to whom correspondence should be addressed.
Water 2023, 15(1), 105; https://doi.org/10.3390/w15010105
Submission received: 7 November 2022 / Revised: 12 December 2022 / Accepted: 20 December 2022 / Published: 28 December 2022
(This article belongs to the Section Hydrogeology)

Abstract

:
In recent decades, the Saïss plain, in the northwest of Morocco, has experienced a noticeable increase in water demand due to a very significant population growth and economic development, as well as the climate change effects. With the aim of reaching optimal and dynamic management of these water resources, it is essential to have comprehensive and reliable information on the state of the aquifer systems in the region. To achieve this, we assessed a geostatistical analysis of groundwater level data, and created a multivariate regression model. Indeed, in this study, a spatiotemporal analysis of groundwater depth based on piezometric measurements of 45 wells was carried out for the period from 2005 to 2020. It compares and evaluates eight geostatistical interpolation methods and solves the problem of data gaps of the piezometric measurement by completing the chronological series of the groundwater level between 2005 and 2020 using the ARIMA model. The results demonstrate that the variation in the groundwater level between 2005 and 2020 indicates that the water table level is decreased in certain areas, but it has improved or remained constant in other areas. These results emphasize an urgent need for a dynamic management for the conservation of groundwater resources in certain areas of the region under this study.

1. Introduction

Groundwater is the main source of water supply for domestic, agricultural, and industrial needs in the Saïss plain in northern Morocco. In recent decades, the water demand has increased due to a very significant population growth and economic development, as well as climate change effects, which have placed remarkable pressure on the water resources. Indeed, climate change, having a significant impact on the availability of water resources, may worsen the country’s situation in terms of the sustainable supply of surface and groundwater. Coupled with growing demand, some regions of Morocco are experiencing water scarcity. Since the early 1980s, the study area has experienced continuous drought; as a result, groundwater exploitation has become more intensive, leading to the decline of groundwater levels and the formation of aquifer depressions [1,2,3].
With the aim of optimal and dynamic management of these water resources, it is essential to have complete and reliable information on the state of the aquifer system of the Saïss plain. Measurements of groundwater levels from observation wells are a source of precise information for analyzing the state of water resources; however, these series are not always continuous in time and space and generally contain gaps for various reasons. Therefore, accurate modeling of groundwater levels in unmonitored locations is necessary for better planning and management.
Understanding, analyzing, and studying the behavior of aquifer systems is essential for making any management decision and for optimum exploitation and rational use of water [4]. The basic tool for this analysis is the regionalized map of piezometric level. It usually serves as a reference to the hydrogeological and environmental studies. It allows, among others, the understanding of the morphology, geometry, and hydrodynamics of the aquifer [5].
Geostatistics is a branch of statistics in which the continuity and spatial behavior of phenomena are incorporated into the interpolation process. It is applied in various fields of Earth sciences, mainly regarding groundwater resources. In addition, the use of geostatistics for the processing of groundwater data extends to several fields of study. It has been implemented by several researchers, particularly for the study of groundwater level variation [6,7,8,9,10,11,12], the study of chemical contamination of groundwater [13,14], and the assessment of groundwater quality [15,16,17,18,19].
In Morocco, geostatistical approaches have been used by many researchers in different fields. For instance, the authors of [20] used geostatistics to interpolate the transmissivity in the Cretaceous basin of Errachidia. It was also used by Lahlou et al. [21] to study spatiotemporal variation of groundwater salinization in the Tadla region. Finally, Rochdane et al. [22] and Nouayti et al. [23] used statistical and geostatistical methods to study the spatial variation of groundwater quality. Unfortunately, geostatistical approaches have not been used for groundwater level monitoring in the Saïss plain.
This work initially aims to test and compare several spatial interpolation methods, namely, deterministic methods (the inverse distance weighting method (IDW), the method of radial networks of the basis function (RBF), and the polynomial local interpolation (PLI)), single-variable probabilistic methods (ordinary kriging (OK), universal kriging (UK), and (EBK)) and multivariate probabilistic methods (ordinary cokriging (OCK) and universal cokriging (UCK)) for the spatiotemporal interpolation of groundwater level in the Saïss plain. All comparisons between the models are assessed and carried out based on statistical accuracy indicators provided by cross-validation and statistics. In the second step, we proceed by solving the problem of missing data from the measurement campaigns of the ABHS piezometric network by reconstructing the time series of the groundwater level between 2005 and 2020 using the ARIMA model. Finally, we investigate and analyze the spatiotemporal evolution of the groundwater resources of the Saïss plain. The scope of this work directly affects the economic decision-making, policies, and strategies for sustainable management of water. Indeed, the choice of the optimal method to interpolate the groundwater level in this region is a starting point for future studies to be conducted.

2. Materials and Methods

2.1. Study Area and Data Collection

The Saïss plain is located in the NW of Morocco with an area of about 2200 km2, 95 km length and 30 km width. It occupies a significant part of the Sebou watershed with the two largest cities (Fez and Meknes) and several centers (Figure 1). The Fez–Meknes basin is part of the south Rifin. The study area delimitation was based on the aquifer system extension and not on the watershed limit.
It is considered a vast asymmetrical syncline in the E–W direction which gradually sinks from south to north and straightens abruptly in contact with pre-Rifin ridges. It is characterized by a low topography decreasing from the south to the north and varying between 250 m and 600 m.
The region is drained by many rivers, including Oued Fez and Oued Mekkès, and their tributaries are the most important. The groundwater system consists of two aquifers (unconfined and confined aquifer). The unconfined aquifer is located in the Plio-Quaternary formation, while the confined aquifer is included in the Liassic limestone and dolomite surmounted by a very thick layer of Miocene marls which can reach 1000 m in thickness. The two aquifers communicate through very deep faults such as the Aïn Taoujdate flexure [24,25]. This study concerns the unconfined aquifer, and all the borehole piezometers are observing only this unconfined aquifer. The region also offers considerable opportunities in terms of thermal mineral water resources, with the springs of Sidi Harazem, Moulay Yacoub, and Ain Allah [26]. The climate of the region is Mediterranean, dry and hot in summer and wet and cool in winter, with significant seasonal temperature differences and an average rainfall of 500 mm in Fez and 600 mm in Meknes. The time series of precipitation, temperature, and potential and actual evapotranspiration reconstructed at the annual and monthly scales are presented in Figure 2 and Figure 3 [27].
The monthly average temperatures are high in summer and generally mild in winter. The average values of potential evapotranspiration (PET) are high. According to Thornthwaite’s formula, the average PET fluctuates between 800 and 900 mm/year. The summer season, almost devoid of rain, records the maximum of the PET, with a maximum in July and August. During the winter, the PET values attenuate considerably, especially in December and January. During the summer, low precipitation, high temperatures, and PET determine a deficit climate balance.
The population of the Fez–Meknes basin is around two million inhabitants and is growing continuously by 3 to 5% per year [28]. Agriculture and crafts are the main economic activities in the Saïss plain.
The piezometric data analyzed in this work were collected from the ABHS. The time series of the groundwater level relate to 39 wells in 2005 and 45 wells in 2020. Other thematic data concerning the aquifer (topography, climatology, geology, etc.) were integrated into a geospatial database managed by a Geographic Information System (GIS) for further processing.

2.2. Methodology

2.2.1. Spatial Variability of the Groundwater Resources

Spatial Interpolation

Our research consists of carrying out a study and a comparative analysis of the different spatial interpolation methods, namely, deterministic methods (IDW, RBF, and LPI), monovariable probabilistic methods (OK, UK, and EBK), and multivariable probabilistic methods (OCK and UCK), in order to choose the spatial interpolation method which seems the most suitable and appropriate for characterizing and mapping the spatial variability of the groundwater level. Then, we impute missing data from the time series of the piezometric network to study and analyze the temporal evolution of the groundwater resources of the Saïss aquifer (Figure 4).
Before starting the spatial interpolation, an exploratory analysis of the spatial data and then a data processing and quality control step were carried out. The latter makes it possible to examine the distribution and the spatial correlation of the variable, verify the normality of the measured datasets, and identify the overall trends of these datasets. The different interpolation methods chosen were then compared. All deterministic and probabilistic spatial interpolation methods can be represented as a weighted sum of the measured data [6,12].
The difference between deterministic and probabilistic (geostatistical) interpolation methods is that deterministic techniques (IDW, RBF, and LPI) do not take into account the spatial correlation between data points and do not provide a measure of the accuracy of estimates, while geostatistical interpolation techniques (kriging) take into consideration the spatial autocorrelation of the measured points and provide indicators of the accuracy of the estimates [6,29]. These probabilistic methods can be applied to variables that indicate a certain homogeneity of characteristics in space (stationary variables) and to nonstationary variables with well-marked tendencies in certain directions [30].
Spatial autocorrelation is examined using a semivariogram, which can be defined as half the variance of the difference between pairs of samples with distance h and whose equation is as follows [31]:
γ h = 1 2 N h i = 1 N h Z x i Z x i + h 2 .
When the experimental semivariogram is defined, a model (theoretical variogram) must be fitted to the points. The modeling of semivariograms is an important step that requires choosing the mathematical function and its parameters which best fit the experimental variogram.
In order to choose the optimal model and compare the interpolation methods, cross-validation was used. This consists of predicting the value of each point in the dataset by removing it and relying on the remaining data. Thus, the difference between the measured value and the estimated value can be calculated.

Criteria for Evaluating the Performance of Interpolation Methods

Three statistical performance evaluation indicators were used to verify the prediction accuracy of the developed models. These are the coefficient of determination (R2), the root-mean-square error (RMSE), and the mean error (ME). These measurements are derived from the formulas given below:
R 2 = i = 1 n P i P ave Q i Q ave 2 i = 1 n P i P ave 2 i = 1 n Q i Q ave 2
RMSE = 1 N i = 1 N   Z ^ x i Z x i 2
ME = 1 N i = 1 N   Z ^ x i Z x i
where:
Pave: Average estimated values;
Qave: Average measured values;
  Z ^ x i : Value to be predicted for a given location;
Z x i : Value observed at the location;
n: Number of values used for estimation.

2.2.2. Temporal Variability of the Groundwater Resources

The aim is to explore the groundwater resources of the aquifer through historical data produced by the ABHS [32]. Predicting these groundwater resources represents a real challenge due to the complexity of hydrological mechanisms [33]. Many numerical models have been developed for this task. Some models are based on physical groundwater modeling. Physical models require adjustment of parameters for each situation [34]. They offer an accurate prediction solution but are difficult to generalize on a large scale. For more systematic prediction needs, the use of time series of groundwater level has been used for years by several authors as an essential tool in water resources planning [35,36,37,38].

Reconstruction of Missing Time Series Data

The collected data in the field may be incomplete due to failures during the acquisition of measurements. This compromises the processing of information. Methods, called imputation, are then applied to consolidate the time series with missing data.
Missing data is a big problem in data science. Knowing the types of missing data can help in choosing the appropriate model to use for the reconstructed ones. There are three types of missing data [39,40]:
Missing data that do not depend on any observable variable and any unobservable parameter are then completely randomly missing data. The gaps in data are then considered to be due to hazards. In this case, the analyses carried out are unbiased.
When missing data are related to the value of an external variable but not to the values of the variable with missing data, we have randomly missing data. This is the most classic case.
If the data are missing for some particular reason, then the data are nonrandom missing.
Most of the methods for imputing missing data are univariate, i.e., the imputation of a variable is performed without using the values of the other variables. These univariate methods are well suited to our case. One of the most popular time series forecasting models is the ARIMA model [41,42,43].
For this study, in most cases, the missing data are alternated by observations in the time series, which allowed their imputation by the model with good precision. According to Osman et al. [37], the percentage of missing data and the mechanism of missing data are the main discriminating criteria for choosing the method to use. ARIMA models are suitable for dealing with multivariate time series. These same criteria for choosing methods were identified in the building sector by Es-Sabar et al. [44], and by Semiromi and Koch [38] in the groundwater sector for processing time series containing up to 40% missing values with acceptable precision.
The mathematical writing of the ARIMA models varies from one author to another. Most of the time, the differences concern the sign of the coefficients. XLSTAT [45] is the most commonly found writing, and is used by most software [46].
If we define by {Xt} a series with mean μ, then, if the series is supposed to follow an ARIMA(p,d,q) (P,D,Q)s model, we can write:
Y t = 1 B d 1 B s D X t μ .
Ø B ϕ B s Y t = θ B Θ B s Z t , Z t 0 , σ 2
with
ϕ z = 1 i = 1 p Ø i z i , ϕ z = 1 i = 1 P ϕ i z i .
θ z = 1 + i = 1 q θ i z i , Θ z = 1 + i = 1 Q Θ i z i .
where:
p is the order of the autoregressive part of the model;
q is the order of the moving average part of the model;
d is the differencing order of the model;
D is the differencing order of the seasonal part of the model;
s is the period of the model (for example, 12, if the data are monthly data, and if one noticed a yearly periodicity in the data);
P is the order of the autoregressive seasonal part of the model;
Q is the order of the moving average seasonal part of the model;
ϕ, Ø, θ, Θ are the polynomial coefficients.
The ARIMA model (autoregressive integrated moving average) is a regular type of time series collected at time intervals. Time series forecasting is the basis of many hydrological analyses and decisions in water resource management. Statistical analysis always pursues two objectives: first, to understand which parameters affected the observation of the time series; and second, the projection of future scenarios based on historical data. An ARIMA (p,q,d) time series supposes static conditions, where d is the number of differences, p is the temporal autoregressive term, and q is the number of moving average [46,47]. The ARIMA model generates data from a time series. In this approach, data generation is based on past and present variables. The main problem of this modeling approach is the determination of the factors p, q, and d. In addition, in simulations using the ARIMA model, all input data must be stationary. In this study, to ensure that all input data were stationary, a process by Yang et al. [48] called differentiation was used.
When applied recursively to each successive observation in the series, each new predicted value is calculated as the weighted average of the current observation and the previously predicted observation; the previously predicted observation was calculated from the observed (previous) value and the predicted value before that (previous) value, and so on.
The objective of the modeling is to determine how many autoregressive parameters (p) and moving averages (q) are needed to obtain an effective model. In practice, the number of parameters, p or q, very rarely exceeds two.
These steps were performed individually for each well. Then, this methodology was validated manually using the observed groundwater level measurements.

3. Results and Discussion

3.1. Data Processing and Quality Control

This essential step aims to identify and verify abnormal samples that risk distorting the calculation of spatial continuity, and apply the interpolation that can be the most representative of the phenomenon studied.
Table 1 gathers the summary statistics of the measurements of the groundwater level of the Saïss aquifer for February 2005 and February 2020. The analysis of the shape parameters, namely, the kurtosis and the skewness, provides information on the distribution of the studied variable. The kurtosis (K = 1.7) being less than 3 highlights a less pointed distribution than the normal law. The skewness (S = 0.22) being greater than 0 indicates an asymmetry in the data distribution. The median (M = 528) and the mean (m = 517) have close values. This proves that the data are substantially close to a normal distribution.
The map of measurement points (Figure 5) served as the basis for this exploratory analysis. It provides information on the continuity and spatial correlation of the variable, and, in addition, makes it possible to identify the existence of trends. In this sense, Figure 6 demonstrates that the studied variable (groundwater level) has a certain spatial continuity (autocorrelation). Indeed, the low values are in the northeast of the study area and the high values in the southwest of it. This observation highlights the existence of a decreasing SW to NE trend.
This remark is confirmed through the trend analysis tool (Figure 6a). Indeed, it uses a three-dimensional reference to represent the measured values, where these are projected onto the XZ and YZ planes. At this stage, the line covering the optimal area from the projected points demonstrates the existence of specific trends [49]. When this line is flat, there is no trend. In our case, Figure 6a, representing the February 2020 groundwater level trend tool, reveals the presence of a drift that can be approximated with a polynomial function. This trend is marked in the N45 direction, which corresponds to the main direction of groundwater flow in this study area.
Moreover, the analysis of the experimental variogram, represented in the form of a variographic surface (Figure 6b), makes it possible to qualify the model as isotropic or anisotropic [48]. The surface represents an isotropic pattern when the variability is the same in all directions. In the present case, Figure 6b presents an anisotropic behavior with large variability in the 45° direction. This is unlike the perpendicular directions, where the variability is constant.
In addition, the variographic cloud tool (Figure 6c) shows the calculated values of the variance γ(h) for all pairs of the samples plotted on the y axis against the distance h separating the two locations, which is represented on the x axis. This tool is used to examine the spatial autocorrelation of the observed phenomenon and to detect anomalous values. Indeed, neighboring points (close to zero on the distance axis) had to be more similar (close to zero on the variance axis). Therefore, close samples with large variances may reveal that the data are inaccurate. The analysis of Figure 6c makes it possible to deduce the spatial correlation and the nonexistence of abnormal values in the sample.
Finally, the analysis of the quality of the data showed the following:
-
The existence of a spatial autocorrelation;
-
A distribution of values approaching the normal distribution;
-
The presence of a drift;
-
Anisotropic behavior;
-
The absence of abnormal values.

3.2. Identification of the Optimal Spatial Interpolation Method

After data analysis and processing, the next step is to identify the optimal interpolation method to generate the most accurate regionalized map. Indeed, the literature offers several interpolation methods adapted to various types of data. These methods make it possible to respond in a specific way to the characteristics of these data.
In this context, the purpose of this section is to demonstrate the absence of a universal rule for choosing the optimal method through data analysis and also to show the importance of comparing interpolation methods. In addition, several authors have directly used universal kriging given the nonstationarity of the variable studied, such as [34,50], contrary to Lahlou et al. [21], who considered the variable studied to be locally stationary. In this case, ordinary kriging was used instead of universal kriging. It is in this sense that this comparative study between different mathematical models was carried out. The methodology followed is illustrated in Figure 4. Several unit tests are necessary in order to identify the optimal method to be used in the interpolation of the groundwater level. Indeed, several parameters and mathematical functions were adjusted and tested for each interpolation method. The methods studied are as follows: deterministic methods (inverse distance weighting (IDW), radial basis functions (RBF), and local polynomial interpolation (LPI)), monovariable probabilistic methods (ordinary kriging (OK), universal kriging (UK), and empirical Bayesian kriging (EBK)), and multivariate probabilistic methods (ordinary cokriging (OCK) and universal cokriging (UCK)).
Multivariate “cokriging” methods take into consideration other variables that are well correlated with the first variable in the interpolation of the studied phenomenon. In this case, the topographic elevation was used as additional data, since it has a good correlation with the groundwater level. Indeed, the two variables are 99% correlated (Figure 7a). In addition, the second advantage of using topographic elevation is that it has good spatial coverage.
Cross-validation played an important role in variographic modeling and in determining the optimal interpolation models. The table below makes it possible to compare the results for each method based on the quality indicators: the root-mean-square error (RMSE), mean error (ME), and the coefficient of determination (R2).
According to the table above, ordinary cokriging (OCK) is the best modeling of reality in this study area with this dataset, since it presents the best quality indicators. Indeed, the OCK has the minimum RMSE and the maximum R2 coefficient. The ranking of the optimal interpolation methods from the most accurate (OCK) to the least accurate (IDW) is given in Table 2.
The assessment of the accuracy of the chosen method is based on the study of the degree of correlation between the measured values and the estimated values. This is performed through cross-validation and statistical indicators. Figure 7b illustrates the linear regression plot of the measured values and those obtained by ordinary cokriging with the optimal model selected in addition to the regression line. The coefficient of determination R2 = 0.987 demonstrates the existence of a strong correlation between the observed values and those predicted.
In addition, Delbari et al. [8] compared mono- and multivariate interpolation methods for estimating the groundwater level using a sample of 260 measurements. The second variable used in its calculation is topographic elevation. In this case, the multivariate methods were more accurate than the monovariate methods.
Figure 8 represents the regionalized maps of the water table depth of the Saïss plain estimated by the optimal methods OCK, OK, and EBK. The analysis of this shows that the flow of groundwater takes place globally along the axis from the southwest to the northeast. Indeed, the high values are found in the southwest (787 m), while the low values (380 m) converge towards the northeast of the aquifer near Fez. In addition, this map provides overall information on the flow velocity. It is high in the NE of the aquifer and moderate in the SW.
The spatial distribution of the uncertainty for the OCK, OK, and EBK interpolations is illustrated in Figure 9. These maps demonstrate the reliability of the predictions in the zones neighboring the samples and make it possible to detect the zones where additional wells will have to be added in order to reach a precise regionalized piezometric map. Indeed, the estimation error varies between 1 and 5 m for the OK, between 2 and 6 m for the OCK, and between 6 and 12 m for the EBK around the measured points, while this error is greater than 10 m in areas devoid of real observations. Therefore, it appears that even the optimal method suitable for the input measurements did not make it possible to interpolate the groundwater level with precision over the entire region.
Finally, the results showed that with the 39 measurements in 2005 and 45 measurements in 2020 used as input, ordinary kriging made it possible to generate the regionalized map of the groundwater level. In the context of the present study, cokriging certainly made it possible to have a good prediction from the point of view of the correlation coefficients R2; nevertheless, ordinary kriging is more precise than cokriging from the point of view of the spatial distribution of uncertainty.
This optimal geostatistical method, in conjunction with mapping tools, made it possible to generate piezometric level maps and standard prediction error maps associated with the kriged one. The maps presented in Figure 9 show that the high interpolation errors are located at the extremities of the plain (mainly in the north, northeast, and south). Therefore, new wells can be implanted in those areas that do not experience excessive exploitation and can be monitored with a reduced temporal frequency.
Therefore, the results of the present comparative study prove, on one hand, the importance of comparing methods to identify the optimal interpolation model. On the other hand, they show the importance of optimizing the piezometric monitoring network of the ABHS in order to have optimal coverage of the Saïss plain. This will allow the smooth running of hydrogeological studies, decision-making, and development, with the ultimate goal of sustainable and integrated groundwater resources management.

3.3. Temporal Groundwater Level Modeling

Time series of groundwater levels often have missing values that need to be considered before using them for hydrogeological analysis, especially for numerical groundwater flow modeling applications. It is more important to ensure sufficient and continuous time series.
Temporal groundwater level simulations were applied using the ARIMA model. An ARIMA model tries to statistically model the time series by considering the average groundwater level as a time series. The applicability of this model was proven in previous research [37,42,51,52].

Reconstruction of Time Series

In this section, we first present the dataset we used, as well as the preprocessing that was carried out. We compiled a dataset corresponding to the period from January 2000 to December 2015 for the piezometers of the region, presenting monthly time series of the groundwater level for which there are more than 10% missing data. Among these piezometers, we can distinguish a group of nine which have more than 10% of missing data. The following table lists the data statistics for each well with a time series with missing data and indicates the number of gaps in their data series (Table 3).
This analysis involved 45 wells that make up the piezometric network for monthly monitoring of the Saïss aquifer. The first step is to complete the time series of the groundwater level over the period from January 2005 to January 2020 using the ARIMA model [35,43]. This period contains several gaps (i.e., 29.7% of missing values). The periods with missing data differ for each piezometer. Figure 10a shows the time series of the groundwater level of the wells with more than 10% of missing data in order to illustrate their imputations.
As indicated in the presentation of the ARIMA model, all data must be stationary in the simulations using this model. The differencing method was applied to make all data stationary. After obtaining a stationary time series, different structures of ARIMA models were obtained through a process of test and error. Finally, ARIMA (1,1,1) was selected for the simulation of the groundwater level in the Saïss plain. The entire statistical calculation was applied by the XLSTAT 2022 software. The simulated groundwater levels by the ARIMA model are shown in Figure 10b.
The applicability of the ARIMA model in groundwater resources simulation, based on the determination coefficient (R2), is illustrated in Figure 11. The results demonstrate that the model produces acceptable performance in simulating groundwater levels.
To assess the consistency between the measured piezometric values and those simulated by the ARIMA model, they were represented graphically. Figure 12 illustrates the observed and the imputation of missing data by the ARIMA model of the groundwater level time series. This figure presents the curves of the observed values and the curves of the simulated values by the ARIMA model.
The ARIMA model simulates a curve of the time series of the groundwater level based on the observed data. This simulation calculated, at the same time, the missing data of this series. This figure shows good agreement between the observed and simulated values. It shows trends and fluctuations in groundwater levels.

3.4. Study of the Variation in the Groundwater Resources between 2005 and 2020

Temporal reconstruction of the chronological series of the groundwater level makes it possible to study its variation. The comparison of the groundwater level maps of 2005 and 2020 provides information on the evolution of the water table depth of the aquifer (Figure 13). The groundwater of the Saïss plain flows from the SW to the NE. At the Aïn Taoujdate flexure, the groundwater flow is diverted towards Fez in the NE and towards Meknes in the NW. The average hydraulic gradient is 5‰. For the comparison of the water table depth between 2005 and 2020, we produced the groundwater level difference map between 2005 and 2020. Figure 13 and Figure 14, respectively, illustrate the map of the wells, according to the value of the groundwater difference, and the interpolated map of the groundwater difference values. These maps make it possible to visualize and analyze the areas that have experienced an improvement or a deterioration in the groundwater resources.
Figure 14 clearly shows that the groundwater level has experienced a significant drawdown in several regions. This phenomenon is caused by population growth, which increased by 48% between 1994 and 2014 [28], and the decline in infiltration. Figure 15 shows a decreasing trend in precipitation between 1960 and 2021 in the study region.
The difference in groundwater level between 2005 and 2020 is classified into three categories (Figure 16):
-
Category 1 corresponds to piezometers in which the groundwater level did not undergo a significant variation. For example, for piezometer 566/21, the groundwater level varies slightly between +/−4 m.
-
Category 2 corresponds to piezometers in which the groundwater level increased by approximately 20 m. These piezometers are mostly located in the irrigated perimeters with surface water. The development of this perimeter made it possible, on the one hand, to reduce the use of groundwater resources for irrigation. On the other hand, it made it possible to have an infiltration of this surface water. At the piezometer 237/15, the water table depth experienced significant fluctuations but with a general upward trend of 25 m.
-
Category 3 corresponds to piezometers whose groundwater levels of the aquifer experienced a decrease of approximately 7 m (Piezometer 3362/15). This drop in groundwater resources is due to excessive water pumping. Moreover, the majority of these piezometers are located in agricultural areas irrigated by groundwater resources. In addition, they are far from potential sources of recharge such as rivers.

4. Conclusions

This study shows how the geostatistics can help to improve spatiotemporal variation of groundwater levels, and evaluates the status of groundwater resources using different geostatistical models. The results prove the inexistence of a universal rule for choosing the optimal interpolation model based on the analysis of the values of the input sample. Comparing methods and assessing their accuracy using cross-validation and statistical indicators before generating the regionalized maps and evaluating their accuracy are the best way. This finding follows the same conclusion of many authors. These tools played a decisive role in the choice and validation of the optimal model.
The identification of the optimal spatial interpolation model is the result of the comparison of different interpolation methods. Kriging made it possible to generate regionalized water table depth maps. The results demonstrate that cokriging produces a good performance in the prediction of groundwater levels (based on R2 and RMSE indicators). However, the ordinary kriging is more accurate than the cokriging from the point of view of the uncertainty of spatial distribution.
In 2005, the Saïss aquifer depth monitoring network managed by ABHS consisted of 39 piezometers. The current network contains 45 piezometers; unfortunately, the spatial coverage is still not optimal. This study will contribute to designing a spatially optimized piezometric monitoring network. In order to implement this new spatially optimal network, the new wells can be implanted first in the extremities of the plain which are devoid of monitoring piezometers and present high interpolation errors. These areas do not experience excessive exploitation and can be monitored with reduced time frequency.
To overcome the problem of time series with missing data, which hinders the smooth running of hydrogeological studies and decision-making, the present study provides a solution based on an ARIMA model to solve this problem. It made it possible to reconstruct the time series of the piezometric monitoring network of the ABHS.
The analysis of the variation in the groundwater resources between 2005 and 2020 shows that, in general, the state of the aquifer decreased in certain areas; on the other hand, it improved or remained constant in other areas. These results will enable the ABHS to target areas suffering from lower water levels through oriented actions to conserve groundwater resources and strengthen the groundwater monitoring network with additional piezometers.

Author Contributions

M.E.G. and H.J.O.: conceptualization, methodology, data collection, data preprocessing, investigation, review of literature, and manuscript preparation—original draft. A.L. and H.R.: supervision, research fund acquisition, validation, reviewing, and formal analysis. All authors were involved in the writing of the manuscript and providing final approval for submission. 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

Processed data are available upon request to the lead author.

Acknowledgments

The authors would like to acknowledge the support from the ABHS (Agence du Bassin Hydraulique de Sebou) for providing data and during the field investigations. We acknowledge the fruitful discussions and contributions of many colleagues and collaborators from the ABHS, FST-Fez, ENSA-Fez, and SAP + D. We acknowledge the reviewers and editors for their valuable suggestions and comments that helped improve the quality of this manuscript.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

  1. Hssaisoune, M.; Bouchaou, L.; Sifeddine, A.; Bouimetarhan, I.; Chehbouni, A. Moroccan Groundwater Resources and Evolution with Global Climate Changes. Geosciences 2020, 10, 81. [Google Scholar] [CrossRef] [Green Version]
  2. Direction de la Recherche et de la Planification de l’Eau (DRPE). Eaux Souterraines au Maroc Comment Concilier Satisfaction des Besoins et Développement Durable des Ressources en Eau Souterraines. Unpublished Report. Ministère délégué auprès du Ministre de l’Energie, des Mines, de l’Eau et de l’Environnement Chargé de l’Eau: Rabat, Morocco, 2014; 1–53. [Google Scholar]
  3. NOVEC-Morocco. Etude pour L’évaluation des Impacts des Changements Climatiques sur les Ressources en Eau et L’identification des Mesures D’adaptation dans le Bassin du Sebou. Unpublished Report. NOVEC-Morocco: Rabat, Morocco, 2020; 200p. [Google Scholar]
  4. Yang, H.F.; Meng, R.F.; Bao, X.L.; Cao, W.G.; Li, Z.Y.; Xu, B.Y. Assessment of water level threshold for groundwater restoration and over-exploitation remediation the Beijing-Tianjin-Hebei Plain. J. Groundw. Sci. Eng. 2022, 10, 113–127. [Google Scholar] [CrossRef]
  5. Chenini, I.; Mammou, A.; El May, M. Groundwater Recharge Zone Mapping Using GIS-Based Multi-Criteria Analysis: A Case Study in Central Tunisia (Maknassy Basin). Water Resour. Manag. 2020, 24, 921–939. [Google Scholar] [CrossRef]
  6. Adhikary, P.P.; Dash, C.J. Comparison of deterministic and stochastic methods to predict spatial variation of groundwater depth. Appl. Water Sci. 2017, 7, 339–348. [Google Scholar] [CrossRef] [Green Version]
  7. Ahmadi, S.; Sedghamiz, A. Geostatistical Analysis of Spatial and Temporal Variations of Groundwater Level. Environ. Monit. Assess. 2007, 129, 277–294. [Google Scholar] [CrossRef]
  8. Delbari, M.; Motlagh, M.B.; Amiri, M. Spatio-temporal variability of groundwater depth in the Eghlid aquifer in southern Iran. Earth Sci. Res. J. 2013, 17, 105–114. [Google Scholar]
  9. Khairul, H.; Sondipon, P.; Tareq, J.C.; Anzhelika, A. Analysis of groundwater table variability and trend using ordinary kriging: The case study of Sylhet, Bangladesh. Appl. Water Sci. 2021, 11, 120. [Google Scholar] [CrossRef]
  10. Khazaz, L.; Oulidi, H.; Moutaki, S.; Ghafiri, A. Comparing and Evaluating Probabilistic and Deterministic Spatial Interpolation Methods for Groundwater Level of Haouz in Morocco. J. Geogr. Inf. Syst. 2015, 7, 631–642. [Google Scholar] [CrossRef] [Green Version]
  11. Manzione, R.L.; Castrignanò, A. A geostatistical approach for multi-source data fusion to predict water table depth. Sci. Total Environ. 2019, 996, 133763. [Google Scholar] [CrossRef]
  12. Xiao, Y.; Gu, X.; Yin, S.; Shao, J.; Cui, Y.; Zhang, Q.; Niu, Y. Geostatistical interpolation model selection based on ArcGIS and spatio-temporal variability analysis of groundwater level in piedmont plains, northwest China. SpringerPlus 2016, 5, 425. [Google Scholar] [CrossRef] [Green Version]
  13. Melloul, A.; Boughriba, M.; Boufaida, M. Étude de la contamination des ressources en eaux souterraines et cartographie de la vulnérabilité d’un aquifère soumis au climat semi-aride méditerranéen: Cas de la plaine côtière de Saïdia, Maroc. Sci. Changements Planétaires Sécheresse 2009, 20, 223–231. [Google Scholar] [CrossRef]
  14. Iskandar, I.; Koike, K. Distinguishing potential sources of arsenic released to groundwater around a fault zone containing a mine site. Environ. Earth Sci. 2011, 63, 595–608. [Google Scholar] [CrossRef]
  15. Nas, B. Geostatistical approach to assessment of spatial distribution of groundwater quality. Pol. J. Environ. Stud. 2009, 18, 1073–1082. [Google Scholar]
  16. Nas, B.; Berktay, A. Groundwater quality mapping in urban groundwater using GIS. Environ. Monit. Assess. 2010, 160, 215–227. [Google Scholar] [CrossRef] [PubMed]
  17. Marko, K.; Al-Amri, N.; Elfeki, A.M. Geostatistical analysis using GIS for mapping groundwater quality: Case study in the recharge area of Wadi Usfan, western Saudi Arabia. Arab. J. Geosci. 2014, 7, 5239–5252. [Google Scholar] [CrossRef]
  18. Mubarak, N.; Hussain, I.; Faisal, M.; Hussain, T.; Shad, M.; Abdel-Salam, N.; Shabbir, J. Spatial Distribution of Sulfate Concentration in Groundwater of South-Punjab, Pakistan. Water Qual. Expo. Health 2015, 7, 503–513. [Google Scholar] [CrossRef]
  19. Venkatramanan, S.; Chung, S.; Ramkumar, T.; Rajesh, R.; Gnanachandrasamy, G. Assessment of groundwater quality using GIS and CCME WQI techniques: A case study of Thiruthuraipoondi city in Cauvery deltaic region, Tamil Nadu, India. Desalination Water Treat. 2015, 57, 12058–12073. [Google Scholar] [CrossRef]
  20. Jarar Oulidi, H.; Benaabidate, L.; Fryar, A.; Benslimane, A. Utilisation du SIG pour la modélisation de la variabilité spatiale des valeurs de la transmissivité. In Incertitude et Envrionnement, fin des Certitudes Scientifiques; Allard, P., Fox, D., Picon, B., Eds.; Edisud: Aix en Provence, France, 2007; pp. 75–86. [Google Scholar]
  21. Lahlou, M.; Ajerame, M.M.; Bogaert, P.; Bousetta, B. Spatiotemporal Variability and Mapping of Groundwater Salinity in Tadla: Geostatistical Approach. In Developments in Soil Salinity Assessment and Reclamation; Shahid, S.A., Abdelfattah, M.A., Taha, F.K., Eds.; Springer: Dordrecht, The Netherlands, 2013; pp. 167–182. [Google Scholar]
  22. Rochdane, S.; Reddy, D.V.; El Mandour, A. Hydrochemical and Isotopic Characterisation of Eastern Haouz Plain Groundwater, Morocco. Environ. Earth Sci. 2015, 73, 3487–3500. [Google Scholar] [CrossRef]
  23. Nouayti, A.; Khattach, D.; Hilali, M.; Nouayti, N.; Arabi, M. Assessment of groundwater quality using statistical techniques in high Basin of Guir (Eastern High Atlas, Morocco). Mater. Today Proc. 2019, 13, 1084–1091. [Google Scholar] [CrossRef]
  24. Essahlaoui, A. Contribution à la Reconnaissance des Formations Aquifères Dans le Bassin de Meknès—Fès (Maroc): Prospection Géoélectrique, Etude Hydrogéologique et Inventaire des Ressources en Eau. Ph.D. Thesis, Mohammadia School of Engineering, Rabat, Morocco, 2000; 250p. [Google Scholar]
  25. Essahlaoui, A. Etude par Prospection Géoélectrique Dans le Plateau de Meknès et Essai de Reconnaissance du Bassin Hydrogeologique du Saïss; Mémoire de DEA, Faculty of Sciences of Meknes: Meknes, Morocco, 1997. [Google Scholar]
  26. Benaabidate, L. Caractérisation du Bassin Versant de Sebou: Hydrogéologie, Qualité des Eaux et Géochimie des Sources Thermales. Ph.D. Thesis, Faculté des Sciences et Techniques de Fès, Fes, Morocco, 2000; 250p. [Google Scholar]
  27. El Garouani, M.; Amyay, M.; Lahrach, A.; Oulidi, H.J. Land Surface Temperature in Response to Land Use/Cover Change Based on Remote Sensing Data and GIS Techniques: Application to Saïss Plain, Morocco. J. Ecol. Eng. 2021, 22, 100–112. [Google Scholar] [CrossRef]
  28. HCP (Haut-Commissariat au Plan, Maroc). Recensement général de la population et de l’habitat au Maroc. 2015. Available online: http://www.hcp.ma (accessed on 1 June 2022).
  29. Antonakos, A.; Lambrakis, N. Spatial interpolation for the distribution of groundwater level in an area of complex geology using widely available GIS tools. Environ. Process. 2021, 8, 993–1026. [Google Scholar] [CrossRef]
  30. Triki, I. Evaluation de techniques d’interpolation spatiale de la piézométrie à l’aide de l’extension Geostatistical Analyst d’ArcGIS. Cas du système aquifère phréatique de Sfax (Tunisie). Géomatique Expert 2014, 99, 55–63. [Google Scholar]
  31. Burgess, T.M.; Webster, R. Optimal interpolation and isarithmic mapping of soil properties I: The semivariogram and punctual kriging. J. Soil Sci. 1980, 31, 315–331. [Google Scholar] [CrossRef]
  32. Agence du Bassin Hydraulique de Sebou (ABHS). Caractérisation de la Nappe de Saïss et Proposition d’un Programme de Renforcement du Réseau de Surveillance des Aquifères; Rapport ABHS, Phase1; Agence du Bassin Hydraulique de Sebou: Fès, Morocco, 2021; 110p. [Google Scholar]
  33. Brédy, J.J.; Gallichand, P.; Gumiere, S.J.C. Water table depth forecasting incranberry fields using two decision-tree-modeling approaches. Agric. Water Manag. 2020, 233, 106090. [Google Scholar] [CrossRef]
  34. Nayak, P.C.; Rao, Y.S.; Sudheer, K. Groundwater level forecasting in a shallow aquifer using artificial neural network approach. Water Resour. Manag. 2006, 20, 77–90. [Google Scholar] [CrossRef]
  35. Beuchée, L.; Guyet, T.; Malinowski, S. Prédiction du niveau de nappes phréatiques: Comparaison d’approches locale, globale et hybride. In Proceedings of the EGC 2022—Conférence francophone sur l’Extraction et la Gestion des Connaissances, Blois, France, 24–28 January 2022; p. hal-03548071. [Google Scholar]
  36. Kisi, O.; Shiri, J.; Nikoofar, B. Forecasting daily lake levels using artificial intelligence approaches. Comput. Geosci. 2012, 41, 169–180. [Google Scholar] [CrossRef]
  37. Osman, M.S.; Abu-Mahfouz, A.M.; Page, P.R. A Survey on Data Imputation Techniques: Water Distribution System as a Use Case. IEEE Access 2018, 6, 63279–63291. [Google Scholar] [CrossRef]
  38. Semiromi, M.T.; Koch, M. Reconstruction of groundwater levels to impute missing values using singular and multichannel spectrum analysis: Application to the Ardabil Plain, Iran. Hydrol. Sci. J. 2019, 64, 1711–1726. [Google Scholar] [CrossRef]
  39. Allison, P.D. Multiple imputation for missing data: A cautionary tale. Sociol. Methods Res. 2000, 3, 301–309. [Google Scholar] [CrossRef] [Green Version]
  40. Benneuil, N. Traitement des données manquantes dans les séries issues des registres paroissiaux. Popul. Hist. 1998, 1–2, 249–270. [Google Scholar] [CrossRef]
  41. Belkacem, B. Prévision avec modèle ARIMA des températures mensuelles du bassin versant de la Seybouse du Nord Est Algérien. Rev. Sci. Tech. LJEE 2018, 32–33. [Google Scholar]
  42. Choubin, B.; Malekian, A. Combined gamma and M-test-based ANN and ARIMA models for groundwater fluctuation forecasting in semiarid regions. Environ. Earth Sci. 2017, 76, 538. [Google Scholar] [CrossRef]
  43. Delcor, L.; Delignières, D.; Cadopi, M.; Brouille, D. Analyse de la variabilité par les modèles ARIMA: Une source d’information pour la compréhension des processus mnésiques. L’Année Psychol. 2008, 108, 699–720. [Google Scholar] [CrossRef]
  44. Es-Sabar, A.; Pannier, M.L.; Godon, A.; Bigaud, D. Traitement des données manquantes pour des capteurs de bâtiments connectés. In Proceedings of the Conférence IBPSA France, Châlons en Champagne, France, 19–20 May 2022; p. hal-03687624. [Google Scholar]
  45. Addinsoft. XLSTAT User Guide. 2022. Available online: https://www.xlstat.com (accessed on 10 June 2022).
  46. Hyndman, R.J.; Athanasopoulos, G. Forecasting: Principles and Practice. OTexts. Available online: https://otexts.com/fpp3/ (accessed on 10 June 2022).
  47. Jeihouni, E.; Mohammadi, M.; Ghazi, B. Response of the Shabestar Plain aquifer to climate-change scenarios through statistical and hybrid soft computing techniques. Groundw. Sustain. Dev. 2021, 15, 100649. [Google Scholar] [CrossRef]
  48. Yang, Q.; Wang, Y.; Zhang, J.; Delgado, J. A comparative study of shallow groundwater level simulation with three time series models in a coastal aquifer of South China. Appl. Water Sci. 2018, 7, 689–698. [Google Scholar] [CrossRef] [Green Version]
  49. Kumar, A.; Maroju, S.; Bhat, A. Application of ArcGIS geostatistical analyst for interpolating environmental data from observations. Environ. Prog. 2007, 26, 220–225. [Google Scholar] [CrossRef]
  50. Kumar, V. Optimal contour mapping of groundwater levels using universal kriging—A case study. Hydrol. Sci. J. 2007, 52, 1038–1050. [Google Scholar] [CrossRef]
  51. Sheikh Khozania, Z.; Banadkookib, F.B.; Ehteramc, M.; Najah Ahmedd, A.; El-Shafiee, A. Combining autoregressive integrated moving average with Long Short-Term Memory neural network and optimisation algorithms for predicting ground water level. J. Clean. Prod. 2022, 348, 131224. [Google Scholar] [CrossRef]
  52. Valipour, M.; Banihabib, M.E.; Behbahani, S.M.R. Comparison of the ARMA, ARIMA, and the autoregressive artificial neural network models in forecasting the monthly inflow of Dez dam reservoir. J. Hydrol. 2013, 476, 433–441. [Google Scholar] [CrossRef]
Figure 1. Geographic location map of the study area.
Figure 1. Geographic location map of the study area.
Water 15 00105 g001
Figure 2. Ombrothermic diagram for the Saïss plain (1988–2021).
Figure 2. Ombrothermic diagram for the Saïss plain (1988–2021).
Water 15 00105 g002
Figure 3. Climate water balance in the Saïss plain (1988–2021).
Figure 3. Climate water balance in the Saïss plain (1988–2021).
Water 15 00105 g003
Figure 4. Methodology flowchart.
Figure 4. Methodology flowchart.
Water 15 00105 g004
Figure 5. Distribution map of used piezometers in 2020.
Figure 5. Distribution map of used piezometers in 2020.
Water 15 00105 g005
Figure 6. Control of data quality. (a) Trend tool. (b) Variographic surface: an anisotropic behavior with variability in the 45° direction. (c) Semivariographic cloud: calculated values of the variance γ(h) for all pairs of the samples plotted on the y axis against the distance h separating the two locations, represented on the x axis.
Figure 6. Control of data quality. (a) Trend tool. (b) Variographic surface: an anisotropic behavior with variability in the 45° direction. (c) Semivariographic cloud: calculated values of the variance γ(h) for all pairs of the samples plotted on the y axis against the distance h separating the two locations, represented on the x axis.
Water 15 00105 g006
Figure 7. (a) Correlation between the groundwater level and the topographic elevation. (b) Correlation between measured and predicted groundwater levels for 2020.
Figure 7. (a) Correlation between the groundwater level and the topographic elevation. (b) Correlation between measured and predicted groundwater levels for 2020.
Water 15 00105 g007
Figure 8. Groundwater level map interpolated by ordinary cokriging (OCK), ordinary kriging (OK), and empirical Bayesian kriging (EBK).
Figure 8. Groundwater level map interpolated by ordinary cokriging (OCK), ordinary kriging (OK), and empirical Bayesian kriging (EBK).
Water 15 00105 g008
Figure 9. Groundwater level prediction errors map for ordinary kriging (OK), ordinary cokriging (OCK), and empirical Bayesian kriging (EBK).
Figure 9. Groundwater level prediction errors map for ordinary kriging (OK), ordinary cokriging (OCK), and empirical Bayesian kriging (EBK).
Water 15 00105 g009
Figure 10. (a) Original groundwater level time series containing missing values. (b) The reconstructed groundwater levels and gap values filled using ARIMA model for nine selected piezometers.
Figure 10. (a) Original groundwater level time series containing missing values. (b) The reconstructed groundwater levels and gap values filled using ARIMA model for nine selected piezometers.
Water 15 00105 g010
Figure 11. Scatterplot of predicted vs. observed values of groundwater level in six selected piezometric stations.
Figure 11. Scatterplot of predicted vs. observed values of groundwater level in six selected piezometric stations.
Water 15 00105 g011
Figure 12. Comparison between measured and reconstructed groundwater level using ARIMA model for the six borehole piezometers.
Figure 12. Comparison between measured and reconstructed groundwater level using ARIMA model for the six borehole piezometers.
Water 15 00105 g012
Figure 13. Map of wells according to the value of the groundwater level difference between 2005 and 2020.
Figure 13. Map of wells according to the value of the groundwater level difference between 2005 and 2020.
Water 15 00105 g013
Figure 14. Interpolated map of the groundwater level difference values between 2005 and 2020.
Figure 14. Interpolated map of the groundwater level difference values between 2005 and 2020.
Water 15 00105 g014
Figure 15. Decreasing trend in precipitation between 1960 and 2021 in the study region.
Figure 15. Decreasing trend in precipitation between 1960 and 2021 in the study region.
Water 15 00105 g015
Figure 16. Hydrographs of the groundwater level for three borehole piezometers.
Figure 16. Hydrographs of the groundwater level for three borehole piezometers.
Water 15 00105 g016
Table 1. Summary statistics of measurements of the groundwater level of the Saïss aquifer for 2005 and 2020.
Table 1. Summary statistics of measurements of the groundwater level of the Saïss aquifer for 2005 and 2020.
DateMinMaxAverageMedianStandard
Deviation
SkewnessKurtosis1st Quartile3rd Quartile
20053737625285171220.221.7406621
20203807875365281190.231.8419627
Table 2. Comparison between optimal models of spatial interpolation.
Table 2. Comparison between optimal models of spatial interpolation.
2020RMSER2ME
OCK13.250.987−0.49
UCK13.280.984−0.48
EBK13.770.9790.19
OK16.040.978−0.86
UK16.140.9770.97
LPI17.710.910−0.94
RBF37.120.900−1.51
IDW45.650.860−9.8
Table 3. Piezometers with gaps in their time series of groundwater level.
Table 3. Piezometers with gaps in their time series of groundwater level.
PiezometerValid MeasurementsMissing Data%
237/151068644.8
787/211147840.6
566/211167639.6
1309/221177539.1
2813/151217137.0
2366/151524020.8
2792/151563618.8
2604/151583417.7
290/221722010.4
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

El Garouani, M.; Radoine, H.; Lahrach, A.; Jarar Oulidi, H. Spatiotemporal Analysis of Groundwater Resources in the Saïss Aquifer, Morocco. Water 2023, 15, 105. https://doi.org/10.3390/w15010105

AMA Style

El Garouani M, Radoine H, Lahrach A, Jarar Oulidi H. Spatiotemporal Analysis of Groundwater Resources in the Saïss Aquifer, Morocco. Water. 2023; 15(1):105. https://doi.org/10.3390/w15010105

Chicago/Turabian Style

El Garouani, Manal, Hassan Radoine, Aberrahim Lahrach, and Hassane Jarar Oulidi. 2023. "Spatiotemporal Analysis of Groundwater Resources in the Saïss Aquifer, Morocco" Water 15, no. 1: 105. https://doi.org/10.3390/w15010105

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop