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

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 km², 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]:

$$\mathsf{\gamma }\left(\mathrm{h}\right)=\frac{1}{2\mathrm{N}\left(\mathrm{h}\right)}{\sum }_{\mathrm{i}=1}^{\mathrm{N}\left(\mathrm{h}\right)}{\left[\mathrm{Z}\left({\mathrm{x}}_{\mathrm{i}}\right)-\mathrm{Z}\left({\mathrm{x}}_{\mathrm{i}}+\mathrm{h}\right)\right]}^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:

$${\mathrm{R}}^2=\frac{{\left[{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{P}}_{\mathrm{i}}-{\mathrm{P}}_{\mathrm{ave}}\right)\left({\mathrm{Q}}_{\mathrm{i}}-{\mathrm{Q}}_{\mathrm{ave}}\right)\right]}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{P}}_{\mathrm{i}}-{\mathrm{P}}_{\mathrm{ave}}\right)}^2{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Q}}_{\mathrm{i}}-{\mathrm{Q}}_{\mathrm{ave}}\right)}^2}$$

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{N}}\sum _{\mathrm{i}=1}^{\mathrm{N}}{\left[ \hat{\mathrm{Z}}\left({\mathrm{x}}_{\mathrm{i}}\right)-\mathrm{Z}\left({\mathrm{x}}_{\mathrm{i}}\right)\right]}^2}$$

$$\mathrm{ME}=\frac{1}{\mathrm{N}}\sum _{\mathrm{i}=1}^{\mathrm{N}}\left[ \hat{\mathrm{Z}}\left({\mathrm{x}}_{\mathrm{i}}\right)-\mathrm{Z}\left({\mathrm{x}}_{\mathrm{i}}\right)\right]$$

where:

Pave: Average estimated values;

Qave: Average measured values;

$\hat{\mathrm{Z}}\left({\mathrm{x}}_{\mathrm{i}}\right)$: Value to be predicted for a given location;

$\mathrm{Z}\left({\mathrm{x}}_{\mathrm{i}}\right)$: 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 {X_t} 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={\left(1-B\right)}^d{\left(1-B^s\right)}^D{X}_t-\mu .$$

$$Ø\left(B\right)\varphi \left(B^s\right)Y_t=\theta \left(B\right)\Theta \left(B^s\right)Z_t,\hspace{1em}{Z}_t\sim \left(0,{\sigma }^2\right)$$

with

$$\varphi \left(z\right)=1-{\sum }_{i=1}^p{Ø}_i{z}^i,\hspace{1em}\varphi \left(z\right)=1-{\sum }_{i=1}^P{\varphi }_i{z}^i.$$

$$\theta \left(z\right)=1+{\sum }_{i=1}^q{\theta }_i{z}^i,\hspace{1em}\Theta \left(z\right)=1+{\sum }_{i=1}^Q{\Theta }_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. Summary statistics of measurements of the groundwater level of the Saïss aquifer for 2005 and 2020.

Date	Min	Max	Average	Median	Standard Deviation	Skewness	Kurtosis	1st Quartile	3rd Quartile
2005	373	762	528	517	122	0.22	1.7	406	621
2020	380	787	536	528	119	0.23	1.8	419	627

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 (R²).

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 R² coefficient. The ranking of the optimal interpolation methods from the most accurate (OCK) to the least accurate (IDW) is given in

Table 2. Comparison between optimal models of spatial interpolation.

2020	RMSE	R2	ME
OCK	13.25	0.987	−0.49
UCK	13.28	0.984	−0.48
EBK	13.77	0.979	0.19
OK	16.04	0.978	−0.86
UK	16.14	0.977	0.97
LPI	17.71	0.910	−0.94
RBF	37.12	0.900	−1.51
IDW	45.65	0.860	−9.8

.

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 R² = 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 R²; 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. Piezometers with gaps in their time series of groundwater level.

Piezometer	Valid Measurements	Missing Data	%
237/15	106	86	44.8
787/21	114	78	40.6
566/21	116	76	39.6
1309/22	117	75	39.1
2813/15	121	71	37.0
2366/15	152	40	20.8
2792/15	156	36	18.8
2604/15	158	34	17.7
290/22	172	20	10.4

).

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 (R²), 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 R² 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.