Homogenization and Trends Analysis of Monthly Precipitation Series in the Fez-Meknes Region, Morocco

: High quality and long-term precipitation data are required to study the variability and trends of rainfall and the impact of climate change. In developing countries like Morocco, the quality of climate data collected from various weather stations faces numerous obstacles. This paper presents methods for collecting, correcting, reconstructing, and homogenizing precipitation series of Morocco’s Fez-Meknes region from 1961 to 2019. Data collected from national specialized agencies based on 83 rain gauge stations was processed through an algorithm specially designed for the homogenization of climatic data (Climatol). We applied the Mann-Kendall test and Sen’s slope estimator to raw and homogenized data to calculate rainfall trend magnitudes and signiﬁcance. The homogenization process allows for the detection of a larger number of stations with statistically signiﬁcant negative trends with 95% and 90% conﬁdence levels, particularly in the mountain ranges, that threatens the main sources of water in the largest watershed in the country. The regionalization of our rain gauge stations is highlighted and compared to previous studies. The monthly and annual means of raw and homogenized data show minor differences over the three main climate zones of the region.

Among the major problems that hinder any study of climate variability are the availability and quality of the data collected by the various measurement stations [11]. The quality of these data in developing countries is often uncertain [12], and their spatial and temporal coverage is limited. The homogenization of climate data is imperative in order to gain a clear vision of climate variability and change, versus changes that are specific to the observation methods [13] (changes in measuring instruments, changes in land use around the stations, changes in the station location, etc.) [14,15]. Climate series are considered to be homogeneous when their variations are related to climate variability alone, and not to these human factors [16]. Changes or alterations in the observation conditions, known as inhomogeneities, can conceal real changes in climate and can mislead the conclusions drawn from their study [17,18].
The detection of these artificial anomalies in climate data is widely studied by scientists. Data on station history (metadata) is one of the methods used to detect the origin of these Largely determined by the relief and exposure, the region of Fez-Meknes has a very varied climate. Different subtypes of Mediterranean climate prevail in the western part, with mean annual rainfall ranging from 423 mm in Fez, 508 mm in Meknes, up to 741 mm in Bab Ounder. The wettest conditions are on the mountain peaks of the Atlas and the Rif (988 mm in Ifrane, 1555 mm in Jbel Oudka). The southeast is arid to semi-arid, especially in the upper Moulouya, with 137 mm at Outat El Haj and 161 mm at Missour. The rainy season starts around October and ends around April, while the dry season lasts between May and September.
The region covers 5.6% of the national territory and is home to 4.2 million inhabitants, roughly 12.5% of the Moroccan population [32]. It provides about 15% of the usable agricultural land in Morocco. The region's contribution to national cereal production remains significant, at around 21% during 2011-2012 [33]. The region's water resources consist mainly of surface water (including the Sebou hydrological basin, which is the richest in water nationally), enhanced by 18 dams, including two large dams (El Wahda dam and Idriss 1st dam). This made it possible to guarantee the irrigation of large agricultural areas, in addition to supplying drinking water to cities not only within the region of Fez-Meknes, but also in the Gharb, which is part of the region of Rabat-Salé-Kénitra through a large project of water transfer [34]. Fossil water is heavily exploited as well. As a result, there are social and environmental imbalances [35][36][37][38]. The increase in water exploitation, combined with the country's recent drought, as well as extreme events such as floods and torrential rains [39], raises the potential for serious problems for the region in the context of climate change [4].

Methods
The homogenization of the precipitation series is based on the work of Guijarro [23] as implemented in the Climatol software tool. The method is based on orthogonal regression, widely known as the reduced major axis (RMA), which is specifically formulated to handle errors in both the x and y variables. It is a popular and alternative strategy to ordinary least squares (OLS) for fitting bivariate relations [40].
The model fitting procedure supposes a linear regression model y = α + βx + ε, where y and x are dependent and independent variables, respectively, α and β are OLS regression coefficients, and ε is a random error term. The RMA slope coefficient is β RMA = β/|r yx |, where r yx is the Pearson correlation coefficient between y and x. The standard error S.E. of β RMA is equal to the S.E. of β. The RMA intercept coefficient is α RMA = y − β RMA , and the S.E. of α RMA is equal to the S.E. of α. [41,42].
This algorithm was utilized through the homogenization process to calculate and estimate the missing data in our series iteratively based on fitting linear relationships between stations in the region until stabilization of the averages is achieved. The algorithm arranges stations in a tree structure, where clusters that show high correlations are grouped in the same branch. The correction of outliers is carried out automatically, and predefined criteria are used to monitor the operation's quality, through the detection of relative inhomogeneities using the Standard Normal Homogeneity Test (SNHT) [19].
As part of Climatol, based on the rainfall behavior of each station, a tree of clusters of stations with similar subregional precipitation behavior was constructed, both before and after the homogenization of data. A generalization map was drawn to present the final regionalization of our stations.
Once the precipitation series have been homogenized and their missing data filled in, precipitation trends were calculated with the Mann-Kendall non-parametric test [43,44], at three levels of significance and the Sen's slope estimator [45]. The goal of this study was to look at the rainfall trend in the region from 1961 to 2019. These results are expected to help in the planning of the water resources management in this area, improving the quality of life of its inhabitants and increasing the resilience in front of the precipitation variability induced by climate change.

Data
The data used in this study was collected from the appropriate administrative authorities in this field, namely the Sebou hydraulic basin agency (ABHS), the Oum Er Rabia hydraulic basin agency (ABHOER), and the hydraulic basin agency of the Moulouya (ABHM), as well as the National Meteorological Directorate (DMN) for some reference stations. Among other things, these agencies look after the water resources in each watershed and gather rainfall statistics. The spatial distribution of the measurement stations is shown in Figure 2, while the other geographic information such as each station's elevation, length of the observation period, and percentage of missing data are shown in Table 1.
The bulk of these stations are within the Fez-Meknes region. There are also some stations just outside this region, to make up for data gaps in places where the rain gauge network is poor within the region. Figure 2 shows the location and code of each station to facilitate the switching between the data table (Table 1) and the map.

Quality Control and Homogenization
Before proceeding with the final homogenization of the series, we used the Climatol package to correct and reconstruct the missing data, where the ratio to the mean is the type of normalization that has been defined with the Climatol tool to allow the comparison of data from stations with very different mean rainfall [23]. After examining the diagnostic graphs of the data generated by Climatol, we increased the threshold to accept outliers up to 14 standard deviations. This was in order to tolerate certain extreme values that represent the region's climatic reality, notably in mountainous areas and during certain seasons when this type of phenomenon occurs (for example, mountain storms during the period of transition between wet spring and summer or between summer and fall).
Referring to the quality flags file generated by the tool (Table 2), the final database consists of 58,764 monthly values for 83 measurement stations between 1961 and 2019, where 61.5% of the data is observed as original and unmodified. Another 4.2% was data that have been corrected after comparison with nearby stations: these are in particular outliers that do not necessarily reflect climate reality. The remaining 34.3% consists of data reconstituted by the Climatol package through adopting the method of [24], a process in which the contribution of nearby stations is critical. The original data in our database is illustrated in Figure 3. The year 1961 has been picked as the start date, to have a beginning relatively well covered with data on the whole of the regional territory, while the end of the data base is dictated by the need to use the latest data available. The amount of information available is more than enough to develop homogenized series and estimate regional and subregional precipitation trends. In Figure 3, the green and red lines respectively mark the desired number (critical threshold) of data availability required for this package. Overall, during the 1960s the amount of data gathered was increasing. After 1970 the network became relatively dense until the 2000s. We then notice the beginning of the observation network's degeneration, which had a negative impact on the quality and quantity of the data collected.
Correlograms illustrate the correlation vs distance (calculated with first differences) in the station series between the raw and the homogenized data. Figure 4 shows a change in favor of homogenized data, where the physically unlikely large negative correlations between a few stations were corrected by the homogenization process. We note that stations less than 150 km apart are better correlated with each other. With increasing distance, a decrease in the correlation coefficient is noticeable in both the original and homogenized data. This prompted us to use the homogenized data to study the final rainfall trends in the region. The detection of breaks in the time series allowed the identification of 13 stations affected by anomalous changes in the mean, including two stations (Tahla and Pont Sebou) with a maximum of two breaks each. These discontinuities are concentrated between the years 1970 and 1980, when Morocco experienced the driest years during the 20th century, so some discontinuities may be due to real changes in climate and not just station-specific factors. This is in agreement with the conclusions of [2,46,47].
The quality of the homogenization applied to the database is also evaluated by the root mean square error (RMSE), between observed and estimated data ( Figure 5). This indicator varies between a minimum of 1.2 mm and a maximum of 58.3 mm, the median is 15.2 mm while the average is 18.3 mm. The standard normal homogeneity test SNHT (on anomaly series) for the adjusted data ranges between a minimum of 0.9 and a maximum of 43.8, the median is 7.4 while the average is 9.0. The comparison between the station's RMSE reveals that the mountain stations register the highest values of the root mean square error. In a mountainous and rugged environment, the challenge of anticipating or interpolating precipitation has previously been noted [48,49].
The comparison of the observed means and those generated after homogenization shows general similarities, and the overall distribution of rainfall remains unchanged ( Figure 6).  Nevertheless, mapping the absolute differences between the means of the two types of data ( Figure 6-bottom) shows a clear change in some mountainous stations, especially the stations with short observation periods, which generate a fairly strong correction through the homogenization process, since the averages observed over a short period do not necessarily represent the local climatic reality (stations Zeggouta S82, Marticha S50, Sidi El Mokhfi S67, Anguid MZ S11, etc.). This correction is also made in special cases of stations located at the top of mountains, such as Tazzeka S78 (1971 m.a.s.l), Jbel Oudka S46 (1589 m.a.s.l), which reflect a very small-scale microclimate.
The comparison of monthly means between the raw series and the homogenized series reflects the strength and validity of the homogenization process. The temporal distribution of the average annual rainfall has also been largely preserved despite the homogenization process.
The absolute difference between the monthly mean values of the two types of data does not exceed the maximum of 4 mm during the month of March, followed by 2.6 mm and 2.3 mm for the months of December and November, respectively (Figure 7). The other values vary between minimums of 0.1 for the months of September and August and 0.4 mm for the months of May and February. The absolute difference between the seasonal averages of the raw data and the homogenized ones remains very small, at 5 mm as a maximum in spring and 1.5 mm in winter. The relationship between the annual averages of the raw data and the homogenized data is shown in (Figure 8).  Across all stations in the region, the maximum annual precipitation occurs in the winter season, with 212 mm during the months (December, January, and February), followed by the spring precipitation with an average height of 156 mm (March, April, and May). The fall season comes in third place with 136 mm. The monthly maximum occurs in the month of December, followed by the month of November, then January, and February. We also notice the increase in precipitation in mountain areas, which is consistent with the findings of Abahous et al. [31].

Regionalization of Stations
Given the station clusters or groups proposed by the Climatol tool (Figure 9), we distinguish three broad zones that are each relatively climatically homogeneous. Their boundaries are largely dictated by the orography, the prevailing wind, and distance to the ocean. The three zones are (1) a humid Atlantic zone with a Mediterranean climate located at the west of the Middle Atlas mountain range; (2) a more arid zone to the south-east and east of the region (Upper Moulouya and Guercif basin); and (3) a transition zone between the two, located on the Atlas peaks, the Medaz basin, and the region's northeastern tip (Figure 9). This transition zone marks the passage between a mild climate shaped by the oceanic influences of the Atlantic through the west of the mountain range of the Middle Atlas and that of the Rif, and an arid climate zone shaped by the combined effects of continentality, the mountain barrier (which induces dry foehn winds), and the great Sahara, which impacts eastern Morocco all the way to the Mediterranean Sea.
A. Humid or "Atlantic" zone This zone is mainly found in the north-western part of the region. Its rainfall regime is presented in Figure 10

B. The Arid zone of Moulouya
This zone is made up of the lands of the upper Moulouya region of Missour, the Guercif basin, and its neighboring areas (Figure 2). The precipitation regime has a very marked spring component with a maximum of rain in April and March, while the second peak happens in autumn, centered on the two months of October and November ( Figure 11). We can distinguish two parts within this zone: The Guercif zone, with more Mediterranean contributions to the north, and the upper Moulouya around Missour, with a double effect on the mountain and the Sahara, which infiltrates via the high plateaus. The absolute difference between the monthly averages of the raw and homogenized data varies between a maximum of 1.7 mm and 1.6 mm during the months of January and November. The minimum is around 0 mm in April and July, and 0.1 to 0.3 mm for the months of February and May, respectively. C. Transition zone or «Mountains» These are the drier mountain peaks of the Middle Atlas and the Rif, which mark the passage between the two previous areas. The variability of precipitation reaches its maximum, and the representativeness of the stations does not necessarily reflect the subregional climatic reality, which is largely shaped by the exposure of the slopes, the opening, and the orientation of the valleys.
The absolute difference between the monthly averages of the raw and homogenized data oscillates between an autumn maximum of 6 mm and 3.6 mm for the months of November and October. The minimum is around 0.3 mm, 0.4, and 1 mm for the months of May, January, and February, respectively, as shown in Figure 12.

Trend Analysis
To examine the trend of precipitation in the region, we have used the Mann-Kendall non-parametric test [43,44] at the alpha thresholds α = 0.01; 0.05; and 0.10, and the Sen slope [45] was applied to the raw and homogenized data.
The results of the trend test are presented in the maps in Figure 13. In both types of data, the trend is negative with significance level α = (0.01; 0.05; 0.10). More of the trends in the raw data are significant at the alpha level α = 0.01 (10 stations). The homogenization process reduces this number to a single station. But this same process made it possible to increase the number of series with negative trends at the SL α = 0.05 from (9 in the raw data to 23 in the homogenized series). For the threshold α = 0.10 the homogenization process increased the number of stations with a negative trend from four in the raw data to 17 stations in the homogenized data. Figure 13. Spatial distribution of significant precipitation trends before (A) and after homogenization (B).
As for the spatial distribution of the negative trends detected, there is a concentration of these negative trends in the mountainous areas of the Middle Atlas and the Rif north of Taza. The arid to semi-arid areas of Moulouya do not show any significant trends, and the same for the western Rif foothills except for two stations (Tabouda and Ratba) which show a downward trend at the 0.10 significance level.
The Sen slope, for the homogenized series which recorded a significant downward trend at the threshold α = 0.05, varies between a maximum of −6.52 and −6, for the Bab Boudir and Echoyeb stations, respectively. Significant minimum values of −2.70, −3.42, −2.38, and −4.05 were seen for the stations of Fez DRH, Meknes, Sidi Kacem, and Taza, respectively. After homogenization, more of the station series thus reflect the global climate trend, which threatens the water resources in the mountains above all, knowing that mountains are the water reservoir for Morocco.

Discussion
In this study, we presented the process of quality control, reconstruction, and homogenization of monthly precipitation data in the Fez-Meknes region, northern Morocco. We also studied the trend of precipitation in the raw and homogenized series between 1961 to 2019.
The comparison between raw and homogenized data through correlograms of the first differences reports a small change related to the correction and reconstruction of the data. The mean monthly rainfall distribution was not greatly affected by homogenization, with differences between the annual and monthly averages of the two data types not exceeding 6 mm.
The quality control flags of the final gap-filled and homogenized dataset showed that our database consists of approximately 62% of the original observed and unmodified data, while 4% is represented by corrected data, and finally, reconstructed data represents nearly 34% of the total.
The validity of the homogenization process is examined through the root mean square error between raw and homogenized data for each station. The RMSE values vary between a minimum of 1.2 mm and a maximum of 80.4 mm, the median is 21.2 mm while the average is 25.3 mm; all these values are far less than those of [31], which may reflect greater station density or more complete data in our study compared to [31]. We have underlined the relevance of the RMSE at stations with shorter observation periods, particularly in mountainous areas where the measuring network is deteriorating. The homogenization operation, therefore, made it possible to correct data from certain stations which represent particular microclimates on the mountain tops but still share temporal variability with nearby stations.
The algorithm's strength is shown in the regionalization of our stations using clusters, which group stations with similar rainfall behavior and divide the Fez-Meknes region into three main climatic zones. These are a humid Atlantic zone to the west, an arid to semi-arid zone to the southeast and east, and finally a transition zone that marks the passage between them. It is worth noting that our results on precipitation regionalization in this Fez-Meknes region are quite similar to those of [50][51][52]; and [47], on the regionalization of precipitation across all of Morocco. There appear to be some differences due to the limited number of stations used by Knippertz [52]. Although the national scale used is much larger than our small region, the three large groups proposed by Ward et al. [50]: the Atlantic zone (I) and the Atlas zone (III), and the Mediterranean zone (V), integrated into the arid zone in our case, are largely recognizable in the regionalization produced in this work.
According to the clustering map ( Figure 9) based on the dendrogram of station clusters, the Middle Atlas stations do not join into a single group, as has been proposed by [47] with its so-called meso grouping. Instead, they are placed in two different groups, one in the Atlantic group for the wettest fringes in the north and west, while the stations located inside the mountain with little rainfall are placed into another group that we have named "Transition zone", towards the arid zones of Moulouya. This corresponds with the so-called micro regions proposed by the same author [47], at least for the Middle Atlas. The northeastern tip of the region also marks the transition between a more humid climate in the west and a drier climate towards the eastern Rif and the Guercif basin in the east.
The power of the Climatol algorithm lies in its ability to detect transition stations between the relative wetness in the west of the region, and the arid zone in the south-east and east. The stations of the Medaz basin, in the eastern and southern parts of Middle Atlas are logically placed in this group that we have called the "transition zone".
The breaks detected by the SNHT test reveal a concentration of ruptures in the 1970s and 1980s, which coincides with a wave of drought which negatively impacted Morocco, in harmony with the results of [2,47,53]. As for the general trends in precipitation series in the Fez-Meknes region, the application of the Mann-Kendall test and the Sen slope to the homogenized series allowed us to observe a general trend towards a decrease in precipitation in mountain areas. This trend is significant at the alpha levels = 0.05 and 0.10. Homogenization made it possible to increase the number of stations with a significant downward trend at the alpha threshold = 0.05, this is in agreement with the results of [7,31,54].

Conclusions
Climate data are important to study climate variability and climate change around the globe, but the quality of this data remains concerning in developing countries. Spatial and temporal discontinuities are the most frequent problems that hinder climate change monitoring. In this study, we presented the first effort at mass collection, reconstruction, and homogenization of rainfall series in the Fez-Meknes region, northern Morocco. Our paper is the second nationally, after that of Abahous et al. [31], to present such an operation in one of the most important regions in the country.
While there are many techniques for homogenization of climatic data, our results suggest that Climatol algorithms perform well in Morocco, consistent with diverse studies around the globe [31,55]. The quality of the result is judged through two elements: the low levels of RMSE between raw and homogenized data, and the regionalization generated for our stations, where data homogenization was accompanied by generating a tree of clusters of homogeneous subregions based on the rainfall behavior of each station. Further, trend analysis in raw and homogenized data shows a clear negative trend in the majority of rain gauge stations located in mountain areas and on the western side of the regional territory. The negative trend of precipitation in mountain areas should prompt decision-makers to rethink adaptation and mitigation strategies in mountain ranges first and foremost. Mountain areas are the main sources of water because of their high levels of rainfall, and many dams and reservoirs are built in those regions. The results of this research will also help in hydrological modeling in the most important watershed of Morocco.