Modeling and Mapping of coastal aquifer vulnerability to seawater intrusion using SEAWAT code and GALDIT index technique: the case of the Rmel aquifer – Larache, Morocco

. The study area of Rmel-O. located in the north of currently faces major water challenges related to the sustainable management of water resources. Climate change and Sea-Level-Rise can increase the risks and costs of water resources management and impact water resources' quantity and quality. Hence, for planning and management, an integrated approach is developed for linking climate models and groundwater models to investigate future impacts of climate change on groundwater resources. Climate projections show an increase in temperature of about 0.45 °C and a reduction in precipitation of 16.7% for 2016-2050. Simulations of seawater intrusion corresponding to various combinations of groundwater extraction predicted climate change and sea-level-rise show that the area will be contaminated on the NW sector of the coastal part. The toe would reach about 5.2 km inland and intrude on high salinity (15–25g/l). Beyond these zones, the contamination of the aquifer will be limited. Moreover, these results were confirmed by the application of the GALDIT method. They reveal that the fringe littoral areas of the aquifer are the most affected by seawater intrusion, with a high risk in the north-western part of the study area.


Introduction
Studying the effect of Climate Change (CC) on groundwater resources is essential to water planners and managers because it may change underlying water management conditions [1,2]. The impact on the groundwater system may be described by referring to its quantity, quality, and functionality. Groundwater recharge is controlled by soil infiltration, land cover, and precipitation amount. Under CC, the variability in the precipitation regime (yearly precipitation amount concentrated in a few precipitation events), combined with land-cover degradation, significantly after groundwater recharge, will be causing a deficit of water balance in the aquifer as well as a degradation of the freshwater quality by seawater intrusion on the coastal plain.
The present study focuses on the Rmel-O. Ogbane coastal aquifer (ROOCA) is intensively exploited. This aquifer is located in the North-West part of Morocco and is very well known for its essential role in the region's industrial, economic, and social development. The paper is based on analyzing the impacts of CC on water resources in the study area based on a 3D groundwater flow modeling and GALDIT method mapping.
The activities can be divided into six main steps: (1) First, analyzing some climate data from available weather stations and collecting future predicted CC datasets with temperature and precipitation variables regarding Rmel-O. Ogbane area in Morocco from 2017 to 2050 based on climate modeling. Indeed, a better understanding of trends can be recognized from observed and forecasted precipitation time series and temperature data in the region, which should be related to CC. More exactly, what impacts these changes will have on the aquifer recharge; (2) based on this climate data (evapotranspiration and rainfall), they are used to estimate recharge for the period 2017-2050; (3) organizing all these data in standard formats stored in a GIS database; (4) creating GIS database better to understand the hydrogeological structure and geometry of the reservoir; (5) identifying and evaluating the spatial distribution of regionalized variables on reservoirs and groundwater resources using geostatistical analysis; (6) analyzing the role of human activities on CC and its impacts on groundwater resources in this area. More specifically, heavy abstraction rates and impacts these changes will have on groundwater level and water quality, and finally (7), both a two-dimensional GALDIT method and a three-dimensional groundwater model are used to identify coastal aquifer vulnerability to seawater intrusion. The developed model considers the increase of groundwater abstractions to satisfy the water demand of urban centers of both Larache and Ksar El Kebir cities. It is capable to better understand the hydrodynamic functioning of the aquifer system by assessing the components of the groundwater mass balance to determine the spatial distribution of the hydrogeological parameters. The model and the method developed will be used as a management tool to assist in planning and water resources management of the aquifer system.
In the upper part of the chart (Fig. 1), tasks assemble several climate data sets for current and future predicted conditions derived from climate modeling, which are used to simulate recharge using the Thornthwaite method. This monthly recharge is then used in a three-dimensional groundwater model (lower part of the chart) to simulate transient saturated variable-density groundwater flow to assess and predict the CC Impacts on groundwater resources and socio-economic vulnerability.
Moreover, the application of the GALDIT method, which consists of mapping the aquifer vulnerability to seawater intrusion (SWI), has been conducted to locate the most affected areas by seawater intrusion.

Background of the study area 2.1 Location of the study area
The Rmel-O. Ogbane aquifer is located in the Low-Loukkos basin immediately at the south of Larache city (Fig. 2). It extends over an area around 305 km 2 and bordures the Atlantic Ocean to the West by the salty level, over a 20 Km band. It is limited by slope gray silts on the left bank of Loukkos wadi and to the East by Mio-Pliocene marl outcrops at the South-East, and the south by the ascended bedrock constitutes a water divide line between the Rmel-O. Ogbane aquifer and the Dradere-Soueire area [3].
The annual average rainfall in the Lower-Loukkos basin reaches more than 1120 mm in the 60 th and declines to 378 mm/y in 2016. The seasonal variation using 12month moving average, the monthly values indicate a decrease over the last five decades after the most important intensity during the 1960s. The variation shows a clear seasonal variability. The monthly values of the temperature in the Lower-Loukkos basin indicate an increase during the last three decades since the 1976's. The impact of CC in the Lower-Loukkos basin cause recurrent droughts, and decreases in recharge directly affect the groundwater level in the aquifer, where it's continuing to decrease. This is associated with excessive over-pumping rates related to the sustainable management of water resources and the delivery of water services for domestic, agricultural, and industrial use. The aquifer supplies good groundwater quality that is easily accessible [4]. This favorable situation has increased pumping and caused environmental problems, such as water table decline and saltwater intrusion on the coast and the coastal plain of the study area. Inappropriate management of this coastal aquifer (overexploitation) may lead to saltwater intrusion into freshwater wells, destroying them as sources of fresh water supply.

GIS Database
A database was developed based on collected data from various studies ( [5,6,7,8]), offices, and organizations such as Direction of the Research and the Planning of Water (DRPE), Hydraulic Basin Agency of Loukkos (ABHL), Regional Office of Agricultural Development of Loukkos (ORMVAL), National Office of Drinking Water of Tangier and Kenitra towns (ONEE), Mohammadia School of Engineers (EMI) and 'Conseil Ingénierie et développement' (CID) which describes the surface water and groundwater resources (hydrology, geology, geometry, groundwater level, and groundwater quality, etc.).
These collected data from the beginning of the 1950s to the year 2016 were collected from the six main organizations mentioned above. All these collected data have been homogenized and stored in a database network of different types: relational geodatabase, geodatabase, and information layers' MONAROO Toolbar' [9].

Hydrogeology and Geostatistics
From the hydrogeological point of view, the Rmel-O. Ogbane aquifer consists of a very heterogeneous complex of fluviolacustrine formations. It consists mainly of Moghrebian shelly sandstones [10] surmounted by Quaternary sands and red silts more or less marly in the Rmel area and laterally evolves in pebbles and sandy silt in the Oulad-Ogbane area. In some places, an intermediate formation of clayey sands and sandy clays is intercalated between the two units of the Moghrebian and Quaternary and can constitute a semi-permeable screen isolating the two hydrogeological units.
A geostatistical model was realized, and some prediction maps were integrated into the GIS. The normality test and the analysis of trends were studied for each variable to select the appropriate Semivariogram and check the results using cross-validation [4].

Conceptual model and boundary conditions
The model grid in the plan view is made up of 128 columns and 128 rows of 250 m (NS-direction) and 250 m (WE-direction) uniform spacing. Moreover, the model grid, vertically oriented, comprises three-layer construction that corresponds to the Rmel hydrostratigraphy layers. Since the model was designed to simulate SWI into both GW flow and solute transfer processes were active in coastal aquifers simultaneously, with boundary conditions taken into consideration. The boundary conditions were set as a constant head boundary around the Atlantic coast for the flow computation, where the cells have a persistent head (mean sea level). The initial head values are set to grid top elevation and turned on in the MODFLOW/SEAWAT program. Moreover, the SEAWAT program was used to simulate the model in a steady-state over 100 years to calculate the initial concentrations of seawater intrusion in 1961. This result will be assigned as the initial concentrations for the transient state. The cells along the Atlantic coast are given a steady concentration of 35 kg/m 3 . Another no line was established through marl outcrops to the southeast of ROOCA and by the risen bottom to the south, which serves as a dividing line between ROOCA and the Dradere-Souiere aquifer. A flow boundary was established in the rest of the aquifer based on the observed data. The aquifer has, for the subsequent layers, internal hydrological stresses are present. The aquifer bottom has a no-flux limit, while the aquifer top has a rebound flood boundary (rain and irrigation return flow).  The simulation, in this case, tests piezometric surface variations over time. Hence, the parameters involved in the temporal equation can be calculated (either a particular yield or specific storage). Aquifer boundary inflows (flux at the boundary, meteoric recharge) and outflows (withdrawals) were also studied. The goal was to replicate the piezometric oscillations found at the control points as accurately as possible during the measurement period. The system reference time is when GW level measurements were made, with the computed piezometric surface in steady-state conditions acting as the initial condition. The model is set up with the initial concentrations of SWI in 1961/1962. The simulation lasted a total of 20075 days, beginning in 1962 and ending in 2016. This simulation makes use of wells that pumping water is mostly used for irrigation and DWS; irrigation started in 1981 and will last until 2016; withdrawals differ over time and are allocated based on piezometric variations. After calibration, the adopted basic storage was estimated to be between 0.1 and 4 percent. This distinction demonstrates adequate consensus between measured and estimated heads in various observation wells for the simulation time under consideration.

Forecasted numerical modeling
The simulation period is based on the groundwater management scenario established by the National Office of Drinking Water (ONEE) in 2016 and ranges from 2017 to 2040. This scenario is based on satisfying future drinking/industrial water demand of both cities Larache and Ksar El Kebir and keeping groundwater abstractions for irrigation constant under CC forcing scenario RCP4.5 and associated with SLR. The numerical model considers three vertical layers to represent the density variation accurately within the aquifer depth. The model was extended one to two kilometers into the Atlantic Ocean, where a constant head and constant salinity boundary was imposed, taking into account the predicted SLR at the case study. Figure 4 shows the results of the predicted drawdown and saltwater extent intruded into the freshwater of the aquifer by 2040.
Under this climate moderate scenario (RCP4.5) and SLR, as integrated into the numerical model, the results show that by 2040, there will be more decrease in hydraulic heads due to intensive pumping discharge to satisfy the water demand of the two urban centers. Indeed, the hydraulic heads will reach negative values around -20 m at the main ONEE well field area, where the drawdown will increase by 25 m, which is highly significant. The aquifer is also contaminated by SWI in the NW sector of the coastal part, in which the toe would reach about 5.2 km inland with an invaded area of about 31 km 2 and will be intruded with high salinity (15-25 g/l) in Layer 3. Consequently, seawater would reach seven observation wells (1534/3, 1535/3, 1536/3, 1380/3, 342/3, 1396/3, and 438/3) and four pumping wells (417/3, 419/3, 718/3, and 1737/3) gradually in 2040. The expected salinity, greater than 2g/l, will be already reached in 2032 and will continue to increase up to 6 g/l in 2040 (Fig. 5). Also, it indicates that the salinity is expected to rise by 2040 from 20 to 30 g/l at piezometer 1380/3, located at 1.5 km far away from the coast.

Results: GALDIT method
The GALDIT method is an intrinsic vulnerability assessment developed by [11] to assess the vulnerability of coastal aquifers to seawater intrusion. It is a parametric system with rating scores and weights for six parameters (Table 1)   The calculation of the GALDIT vulnerability Index is given by the sum of the six index maps at each cell (size: 250 m for our case study) following equation (1): Where R is the rating and W is the weight string, as indicated in table 1. The final GALDIT vulnerability index score varies from 5 to 10 and is divided into three vulnerability classes: high (>10), moderate (7.5 to 10), and low vulnerability (5 to 7.5). The results show that for 2014 (Fig. 6), the vulnerability classes obtained by this combination vary between 5 and 10. They are distributed in three classes corresponding to the degrees of vulnerability fluctuating from "Low" to "High". The dominant class is the "Low vulnerability", covering up to 93% of the study area. From Medium to High Vulnerability classes cover an area up to 5% to 2% respectively and are located in the northwestern part of the study area. Moreover, these results were confirmed by the Numerical modeling of SWI. They reveal that the fringe littoral areas of the aquifer are the most affected by seawater intrusion, with a high risk in the north-western part of the study area.

Conclusion
Groundwater is the important source of freshwater in Rmel O. Ogbane coastal plain in the Low Loukkos basin (Morocco), where 98% of the domestic water and the whole industrial water depend on groundwater. Therefore, it is imperative to protect the groundwater aquifer from seawater intrusion in this area. To predict the extent of seawater intrusion in the future and provide helpful information for the protection of groundwater resources, a GALDIT vulnerability map and a 3D numerical model of density-dependent groundwater flow and miscible salt transport have been developed to assess the current extent of seawater intrusion in the study area. These two approaches are used to better understand coastal aquifer vulnerability to seawater intrusion for the Rmel-O. Ogbane aquifer in this semi-arid region of Morocco under CC and human activities.
The results show that the main effect of seawater intrusion in the Rmel-O. Ogbane aquifer will be the excessive over-pumping that will decrease the renewable resources, which with the pessimist scenario may be increased by 30% of present-day values. Also, the decline in recharge and the rise in sea level due to the CC accelerate saltwater intrusion into the aquifer system, which reduces the fresh groundwater resources. The simulation results show that the maximum extent (about 5.2 km) of seawater intrusion will increase in 2040 in the north-western sector. The water quality will be affected mainly in the ONEE pumping area, immediately adjacent to the seashore.
Moreover, these results were confirmed by the application of GALDIT method. It is capable of assessing the vulnerability to seawater intrusion of the coastal aquifer of Rmel. In 2014, the groundwater aquifer was characterized by a low vulnerability upstream and a high vulnerability downstream with an intense vulnerability of the marine intrusion in the coastal area and the Loukkos River's proximity. According to GALDIT method, the aquifer is highly vulnerable up to 1 km inland in the northwestern part of the study area.
Indeed, this derived 2D and 3D maps can be used as a tool for a decision support system to assist better strategies for groundwater development, planning, and management of the coastal groundwater resources.