Numerical Modeling of Saltwater Intrusion in the Rmel-Oulad Ogbane Coastal Aquifer (Larache, Morocco) in the Climate Change and Sea-Level Rise Context (2040)

: Many coastal aquifers have experienced seawater intrusion (SWI) into fresh groundwater aquifers. The principal causes of SWI include over-pumping and events such as climate change (CC) and rising sea levels. In northern Morocco, the Rmel-Oulad Ogbane coastal aquifer (ROOCA) supplies high-quality groundwater for drinking water and agriculture. This favorable situation has led to increased pumping, resulting in environmental challenges such as dropping water table and SWI. Furthermore, the climate has resulted in less recharge, with an estimated annual precipitation of 602 mm and an average temperature of 18.5 ◦ C. The goal of this study is to determine how CC, over-pumping, and sea-level rise (SLR) affect SWI. Computational groundwater and solute transport models are used to simulate the spatial and temporal evolution of hydraulic heads and groundwater solute concentrations. The calibration is based on steady and transient groundwater levels from 1962 to 2040. SWI simulations show that the NW sector of the coastal area would be polluted, with the toe reaching 5.2 km inland with a signiﬁcant salinity (15–25 g/L). To protect the fresh water in the reservoir from SWI, enhanced groundwater development and management approaches for this aquifer are required, such as artiﬁcial recharge from surface water.


Introduction
Seawater intrusion (SWI) is a worldwide problem that has been exacerbated by rising sea levels, climate change (CC), and excessive over-pumping (EOP) of coastal fresh groundwater (GW) resources. Many of the world's coastal areas are characterized by dense settlements, with a large part of the world's population residing within 60 km from the coast [1]. More than 3 million Moroccans live in Morocco's coastal areas, and this figure is gradually increasing [2]. Indeed, in 2015, more than half of the population lived on the coast, which experienced drought and agricultural development with an increasing proportion of the rural population due to rural flight [2]. As a result, GW overexploitation soon became a widespread problem, with several coastal regions around the world (Morocco [3][4][5][6], Tunisia [7][8][9], Libya [10], Italy [11], and the Netherlands [1]) experiencing substantial SWI in aquifers, resulting in significant degradation of GW quality and, consequently, quantity [12][13][14][15][16].
To better understand the saltwater intrusion process, numerous studies have been undertaken, including laboratory-scale investigations and numerical and analytical modeling studies [17][18][19]. A probabilistic numerical model was built to predict the extent of saltwater intrusion into a coastal phreatic aquifer by Felisa et al. [17]. Fahs et al. [18] used the semianalytical solution of the dispersive Henry problem based on the Fourier technique to simulate heterogeneous and anisotropic coastal aquifers, taking into account the effects of the aquifer hydraulic parameters, boundary conditions, pumping and recharge rates, Ocean warming and ice melting, according to the Intergovernmental Panel on Climate Change (IPCC), could raise global sea levels by up to~60 cm by 2100 [46].
The effects of climate-induced SLR will become more apparent as the intensity of the phenomenon increases, especially in low-elevation coastal areas ( Figure 2). The majority of African countries tend to be in grave danger due to low levels of development along with projections of fast population increases in coastal areas. Under these conditions, the SLR at the aquifer's seaside boundary, which is an additional pressure head, would be imposed. Water table gradient and/or the piezometric head would decrease, resulting in further interference. Many hydraulic, geometric, and transport parameters influence the SWI operation. Each aquifer has its set of circumstances, and the sharp interface approach cannot be applied to all of them, because of the transient conditions which often reign in the exploitation of aquifers. Coastal areas vulnerable to climate-induced sea-level rise. Morocco is indicated in the northwestern part of Africa by the small red frame [47].
For this purpose, the study's goals are to (i) describe and characterize the case study (e.g., geography and climate) (Section 2.1); (ii) define the geological and hydrogeological setting of the study area (Sections 2.2 and 2.3); (iii) understand the hydrodynamic functioning of this hydrogeological system by quantifying the components of the GW mass balance (Section 2.4); (iv) design a conceptual model of the ROOCA (model set-up, initial and boundary conditions, model calibration and validation, time step of the simulations, the distribution of recharge from all sources and pumping, etc.); (v) design a computational model for simulating saltwater intrusion in the ROOCA using SEAWAT code; and (vi) determine the expected future forcing on the ROOCA aquifer due to population increase, water demand (over-pumping), CC (rainfall reduction), and SLR extrapolations to assess the water balance terms and SWI volumes under CC scenarios (Section 3).

Description of the Site
The study site (ROOCA) can be considered representative of the majority of northern Moroccan seashores. It is located in the northwestern part of Morocco in the low Loukkos basin immediately south of Larache city, with an area of about 305 km 2 ( Figure 3). This region is bounded to the west by the Atlantic Ocean along a 20 km band, a succession of hills of Prerifan geological formations to the east, Mio-Pliocene marl outcrops to the southeast, and a risen bottom that acts as a water divide line between the ROOCA and the Dradere-Soueire aquifers to the south [48]. Figure 4 depicts that ROOCA areas are home to more than 4.5 million people. The observed increase in the urban population is attributable to natural demographic growth and the rural exodus from rural to urban space, as well as the formation of new urban centers and the expansion of city borders.  The Atlantic Ocean influences the region's climate, which is subhumid. The monthly temperature values in the Loukkos basin show an increase over the last three decades, dating back to 1976 (Figure 5a). The annual average temperate in July to January ranges from 25.05 to 12.18 • C. State annual rainfall recorded from 1961 to 2016 is around 684 mm. Figure 5b shows a clear seasonal irregularity, where the ombrothermic diagram shows that the majority of the rainfall falls from October to March, with the rest of the year being almost completely dry (Figure 5d). In the study region, the yearly average evapotranspiration is calculated to be 384 mm/year. An influence of CC in the Loukkos basin triggers recurring droughts and reductions in recharge (Figure 5c). This is coupled with EOP ( Figure 5e) to supply industrial and domestic irrigation and agricultural and metropolitan regions. The whole state has experienced a significant decrease in GW level (Figure 5f), which could potentially result in an aquifer water balance deficit and a loss of GW production due to SWI on the study area's coast and coastal plain.    at Larache station, located in the Loukkos basin; (b) precipitation  at Larache station (Loukkos basin); and (c) natural recharge of ROOCA calculated by Thornthwaite method [49]. (d) Ombrothermic chart at Larache station  by Bagnouls and Gaussen [50]. (e) Annual GW pumping volume of ROOCA in hm 3 /year. (f) Decrease in water table derived from the monitoring of two representative wells located in the center (1407/3) and the coast (342/3) of the ROOCA.

Geological Background
The existing geologic units in the research field have been intensively studied [50][51][52] ( Figure 6a). The geological framework of ROOCA is characterized by the Pliocene-Quaternary superposition above a regional, Mio-Pliocene, predominantly marly, layer. It is formed from the bottom to the top by the following layers: (1) The Mio-Pliocene sediments are marked by blue marls, which form the impermeable bedrock of the aquifer.
(2) The Plio-Villafranchian sediments are characterized by coastal and dune deposits that are usually 20 to 50 m thick and are composed of shelly sandstones, sands, and sandy marls. The continental deposits of the Villafranchian period are composed of red-clay cement pebbles surrounded by red sandy clays. (3) The Quaternary: Sandstones and coquina are marine Quaternary remains (ancient Quaternary). Rhamna sands and dune sandstones, as well as numerous fluviatile alluvial deposits, comprise the continental Quaternary and sandstones. The most recent Quaternary alluvial deposits (which contain some blue marls that are more or less sandy and contain marine shells, attributed to the Flandrian transgression) contain primarily red clays with sandy or stony passages, attributed to the Soltanian, and marl-silty layers that are black or greyish, attributed to the Rharbian.
This zone was very sensitive to Post-Villafranchian tectonic movements; they isolated the low Loukkos basin to the north and the Dradere-Souere basin to the south. The structure is characterized by NW-SE orientated ripples showing basins (Rhamna) separated by anticlines (Figure 6a).

Hydrogeology
For the hydrogeological background, Rmel coastal aquifers are made up of Plio-Quaternary sands and sandstone, but the bottom is composed of blue marls [52]. ROOCA is made up of Moghrebian [53] shelly sandstones that are topped in the Rmel field by Quaternary sands and marly red silts. In the Oulad-Ogbane region, they grow sideways in pebbles and coarse silt. Its thickness varies between 0.1 and 146 m. ROOCA's thickest coat is found in a large bowl of Rhamna (140 m), which gradually reduces closer to the edges. The isopach map of the intermediate layer, which comprises clayey soil and sandy clay in the Rmel zone, forms a semipermeable screen that hydraulically isolates the two aquifers (upper and lower) [54]. The ROOCA aquifer is partly confined (limited) and partly unconfined for the rest of the area. These field observations have been taken into account when modeling the aquifer; they are based on pumping test measurements (DRPE, 1987) and show (1) existence of a limited "semiconfined" area with a constant delayed flow (leakage) and (2) existence of an unconfined area with a nonconstant delayed flow. Three geological cross-sections are given and illustrate more details of the existing hydrogeological units: were packed with alluvial deposits and aeolian sands. Shell sandstones resurface in the boreholes 470/3, 700/3, and 704/3 but are missing in borehole 705/3, where the Villafranchian pebbles are found above the Pliocene marl sands [51].
In this analysis, we attempted to model the subsurface of this region in a geoscientific information system (GSIS) using the RockWorks program and borehole data. This program was designed for data collection, viewing, and analysis and displays the reservoir's hydrogeological composition, extension, and geometry ( Figure 7). This software platform is used for incorporating 3D geological models of sedimentary media into traditional hydrogeological modeling methods.
The 127 complete boreholes georeferenced in GIS were integrated into the RockWorks 'Borehole Manager' module that contains borehole processes such as data entry, management, and review ( Figure 7a). The location of these boreholes is illustrated in Figure 7b. The geological units encountered are explicitly mentioned in the geometric model, as shown in Figure 7b-d below. The 'Stratigraphy model' module of RockWorks was used to convert all boreholes data and cross-sections into a 3D model representing the geometric model as a whole. The program creates a grid pattern of each surface and its stratigraphic base. The 3D solid allowed us to determine the geometry, volume, and location of lithostratigraphic deposits in the GW environment of the ROOCA (Figure 7c), and the fence diagram visualizes the configuration between overlapping surfaces (Figure 7d). Additionally, a geostatistical analysis was performed, and some prediction maps were incorporated into the three-dimensional geoscientific information system (3D GSIS). For each measure, the normality test and pattern analysis were used to pick the relevant semivariogram and cross-validate the results [55].
This 3D modeling allowed the design of a conceptual model of the aquifer reservoir, the performance of a set of simulations of GW flow in steady and transient states, and the design of a solute transport model for SWI in the ROOCA.

Water Balance
The regional GW flow is primarily SW-NE, with discharge to the Atlantic Ocean. However, as a result of GW during pumping, a cone of depression shifts the flow pattern to the north. Variations caused by surface water recharge from irrigation and pumping wells can be seen in piezometric data records from 1980 to 2016. The largest renewable groundwater recharge source comes from rainfall and drainage return discharge. Based on a field investigation conducted in 2013/14, the overall GW withdrawal from the pumping wells is evaluated to be 479 L/s; recharge to the aquifer is predicted to be 1508 L/s (Table 1). In 2013/14, the expected outflow to the Atlantic Ocean is 293 L/s. According to the results of a hydrogeochemical analysis conducted between 1985 and 1992, GW salinity in certain coastal observation wells, where measurements were made between the water table level and the aquifer bottom, ranges from 1 to 6 g/L.  The development of a GIS database [56,57] helped to improve and update the ROOCA's hydrogeological water balance based on a mass-balance model with respect to 1963 and 2014 (Table 1 and Figure 8). This was based on a global and comprehensive inventory of all pumping points (wells, boreholes, and springs) and rain data that are reasonably reliable to estimate lateral flow in 2013/14. Natural GW recharge was calculated using the disparity between inflows (rainfall, irrigation) and outflows (plant evapotranspiration, surface run-off, and drinking water supply (DWS)) in a soil water balance equation, as follows:
Assessment of GW withdrawals for rural, DWS, and industrial purposes using indirect calculation or data collection; 3.
Calculation of aquifer inflow and outflow values for nearby water sources.

Model Discretization
Using the database developed by [56,57], the mathematical model developed simulates transient GW flow and solute transport for variable density from 1962 to 2040.
The model grid in the plan view is made up of 128 columns and 128 rows of 250 m (NS direction) and 250 m (S direction) with uniform spacing (WE direction), in that order ( Figure 9a). The model grid, vertically oriented, is made up of three-layer construction that corresponds to the above-mentioned ROOCA hydro-stratigraphy layers (Figure 9b,c). The cross-section through the Rmel coastal plain, as well as the hydrogeologic units and a vertical view of spatial discretization for the quasi-three-dimensional numerical model, is shown by the red line. An inactive area is represented by the green color, while the active area is represented by the white color. The Quaternary aquifers are represented by Layer A1. Layer A2a is an established aquitard between layers A1 and A3 of the Rmel aquifer, and it is thought to be semiimpermeable. When the thickness of Layer A2a exceeds 30 m, Layer A3 reflects the Moghrebian aquifer, which acts as a semiconfining to confining layer. In the Oulad-Ogbane sector, layer A2b represents pebbles and sandy silt. The bottom boundary (A4) in this model is set to Mio-Pliocene. The thickness of the model layer is calculated by (i) the presence of sand and shelly sandstone deposits documented in well driller's documents archived at the Loukkos Hydraulic Basin Agency (ABHL) and the Direction of Water Research and Planning (DRPE) and (ii) a 3D GSIS (Figure 7), which depicts the reservoir's hydrogeological composition, expansion, and geometry. The model simulation used a total of seven stress periods, ranging from 1961/1962 to 2039/2040, with a one-year time step.

Parameters in Hydrogeology
The hydraulic conductivity distribution was determined using data from pumping tests conducted at the time of good completion and recorded in well driller logs [54]. Each zone's hydraulic conductivities are standardized. Pumping tests yielded hydraulic conductivities ranging from 9 × 10 −5 to 2.2 × 10 −4 m/s. For model calibration, the hydraulic conductivities of the model were modified during model tuning to account for the scale of the field. The ABHL administration provided monthly precipitation data. Using the Thornthwaite method, we determined an average evapotranspiration rate of 383.7 mm/year. As a starting point, effective porosity of between 0.15 and 0.25 was chosen. The dispersion coefficients were calculated using a widely used approach based on the sample area's scale [58]. The longitudinal dispersivity was set to 10 m. The ratio of horizontal transverse dispersivity to longitudinal dispersivity (α T /α L ) was thought to be about 0.01. However, it was believed that the vertical transverse dispersivity to longitudinal dispersivity ratio (α V /α L ) was 0.001. The sample region's aquifer property values were assigned to the model cells ( Table 2). The calibration of the model was completed to achieve findings that were as close to the observed concentration data as possible. Table 2. Aquifer parameters used for the study area.

Initial and Boundary Conditions
Since the model was designed to simulate SWI into aquifers, both GW flow and solute transport processes were coupled at the same time, with initial and boundary conditions taken into consideration. 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 steady state over a 100-year period in order to calculate the initial concentrations of SWI in 1961. This result will be assigned as the initial concentrations for the transient state. For the flow computation, the boundary conditions ( Figure 9a) were set as a constant head boundary with a blue line. Around the Atlantic coast, the cells have a persistent head (mean sea level). A no-flow boundary was established through marl outcrops to the southeast of ROOCA and by the risen bottom to the south, which serves as a dividing flow line between ROOCA and the Dradere-Souiere aquifer. A flow boundary was established in the rest of the aquifer (orange color) based on the observed data. The drain boundary is defined as a green line. The aquifer has internal hydrological stresses which are applied for the subsequent layers. The aquifer bottom has a no-flux boundary, while the aquifer top has a recharge boundary (rain and irrigation return flow). The cells along the Atlantic coast are given a specific concentration of 35 kg/m 3 .

Climate and Sea Level Data
Datasets used for this assessment comprise a combination of regional climate modeling projection data generated from Regional Initiative for Assessing Climate Change Impacts on Water Resources and Socioeconomic Vulnerability in the Arab Region (RICCAR) [59] and a set of local observation datasets for precipitation and temperature relevant to our study area. This section is based on extracting time series of Pr (Precipitation) and Ta (Temperature) variables for the entire time period from 1951 to 2100. The available files in the NetCDF format (.nc) were used to extract these time series for our study area located at latitude = 35.2 • North and longitude = 6.16 • West.
Based on RICCAR data, we could present some plots of time series that summarize and show the updated knowledge on the climatology of the study area. We have extracted climate data and plotted some diagrams in Figures 10 and 11, which show the evolution of projected P and T for various climate models and scenarios. The main trends of the parameter variation are also provided to analyze and measure the CC tendency of these parameters. Temperature and precipitation projections are derived from CNRM-CM5, EC-EARTH, and GFDL-ESM2M regional climate models (RCMs) under Representative Concentration Pathway (RCP) 4.5 and RCP 8.5 scenarios. In all cases, the forecasts for these three RCM models indicate a drop in precipitation ( Figure 10) and a rise in temperature ( Figure 11). The main trend for RCP 8.5 is relatively much stronger, as temperature increases more and precipitation decreases more. Hence, CC surely will harm the aquifer recharge and SLR, which may pollute fresh GW coastal areas and reduce their potential recharge.  Based on National Oceanic and Atmospheric Administration (NOAA) data, the following figure (Figure 12) is plotted and provides time series for estimates of SLR based on measurements from satellite radar altimeters such as TOPEX/Poseidon (T/P), Jason-1, Jason-2, and Jason-3, which have been in use since 1991. Glacial isostatic change impacts on the geoid, which are estimated to be +0.2 to +0.5 mm/year when globally combined, are not included in the SLR calculations. For our study area, we then deduce that the sea level will rise from~29.71 mm by 2010 to~71.72 mm by 2020 ( Figure 13). The impacts of climate-induced SLR will become more evident as the severity of the phenomenon grows, especially in low-elevation coastal zones. Besides, the impact of CC (projected temperature and precipitation obtained from RCMs under RCP 4.5 and RCP 8.5 from 2020 to 2100) will have a negative impact on SLR in the study as a result of the melting of glaciers and the warming of the oceans. For this purpose, linear regressions were used to project SLR (yearly values) from 2020 to 2040 from the observed data in answer to the melting of glaciers and the warming of the oceans. This modification of SLR is implemented in our modeling step of GW flow and SWI into the aquifer.

Numerical Model
The GW numerical flow design was created based on using the finite difference method in Visual MODFLOW Flex and includes MODFLOW code [60,61]. In two-and three-dimensional problems, modular finite difference is used in (3D) modeling to solve the differential equation that regulates flow in a porous medium using the GW structure in steady and transient states of the ROOCA.
A 3D numerical GW flow model, calibrated under steady and transient states resolving Equation (1) where K xx , K yy , and K zz are the K-values along of the x, y, and z align directions (L/T); h is the hydraulic head (L); W is the volumetric flux per unit volume and represents sources and (or) sinks (T −1 ); S s is the specific storage of the porous medium (L −1 ); and t is the time.
After that, a 3D variable-density GW flow model, calibrated under steady and transient states resolving Equation (2) for the years 1962 to 2040, was generated using Visual MODFLOW and the SEAWAT code [62,63], which is a valuable method for simulating different variable-density fluids moving through dynamic-geometry hydrogeological environments, such as SWI in coastal aquifers. SEAWAT was created by integrating an updated version of MODFLOW with MT3DMS in a single programming application.
As flow and transport are inextricably linked, they were used to study SWI in the ROOCA under a distribution of salt concentration values in the coastal area.

∂C ∂t
Here, C is the contaminant concentration in GW (ML −3 ), x i is the position in the cartesian coordinate axes (L), D ij is the hydrodynamic dispersion coefficient (L 2 T −1 ), V i is the fluid velocity (LT −1 ), q s is the flow per unit of injected aquifer length (or pumped), Cs is the concentration of recharge or discharge flow (q s ) (ML −3 ), θ is the porosity of the porous medium, and Rk (k = 1, . . . ,N) is the rate of solute production or decay in reaction k of N different reactions (M L −3 T −1 ).

Calibration and Model Results
The main aim of the calibration stage is to produce results that are as close to the field data as possible by modifying the system's parameters. Calibration of the model was accomplished by changing the distribution and values of two main parameters, namely the hydraulic conductivity and the specific storage coefficient of the aquifers. The hydraulic heads were determined by the model until they reach a suitable accuracy to match the observed values.
At first, the model was used in steady-state conditions for its calibration using observed piezometric data from 1961/62 [51]. The few known hydraulic property values were used as input parameters, and simulations were run by adjusting the hydraulic conductivity to obtain the best match between predicted and calculated piezometric values at the available control observation wells.
Following that, the PEST software was used to configure the model's parameters and obtain optimal calibrations for various starting conditions. As seen in Figure 14a,b, the results indicate satisfactory agreement between computed and observed heads. Figure 14b shows the good match between measurements and calculations in the 24 observation wells within the modeled region; the model calibration yielded two sets of values that are highly correlated, as illustrated by the map, which is very similar to the perfect correspondence axis in the scatter diagram, with an average correlation coefficient of 0.999, a mean error of 0.19 m, and a root mean square error of 0.849 m. The mean percentage difference is about 0.827 m, which is a reasonable value given that the computational model cannot account for local differences in the real world, and hydraulic head values were measured using altimetry data collected on a 1:50,000 scale (regional technical cartography). However, this means that the measured and computed heads are relatively similar. The overall model is then closer to reality, with respect to gradients and heads that are almost identical. The residuals between the observed and estimated heads are also indicated in Table 3. The lateral hydraulic conductivity of layers A1, A2, and A3 as a result after model calibration varied from 4.10 −4 to 8.10 −7 m/s (Figure 15b).
The cross-sections and 3D representation of the stratigraphic model indicated the various aquifer levels (hydrogeological units) based on three geological ages ranging from the Plio-Quaternary to the Moghrebian aquifer. However, the lithostratigraphic formations are heterogeneous, varying in each unit where sands, marls, and silts (e.g., upper unit) can be found, which is why the Plio-Quaternary aquifer formation has 11 categories of permeabilities (Table 4). Hence, the hydraulic conductivity varies to account for the lithostratigraphy of the geological formations.
In the ROOCA aquifer, a geostatistical approach was used to study the spatial distribution of regionalized variables based on the pumping test data [55]. This K-value distribution was taken as initial values assigned to the model (initial permeability) and could not produce a simulated piezometry comparable to the measured one (1961/62). To overcome this situation, we used the K-value distribution based on earlier research conducted by Larabi on the same case study [64].      Table 5 shows the observed and calculated water balance with its various elements, along with the release to shore, which is found to be 584 L/s. The key input variable is precipitation recharge, which represents 96% of inflows, and the main output component is drainage to rivers and towards the neighboring alluvial plain, which represents 56% of outflows.
However, in a transient state, the model was calibrated until the end of 2009 due to a lack of consistently controlled hydraulic head performance using observed hydraulic heads from 1961/62. In total, 22 observed hydraulic head simulations were performed during calibration. Figure 16 depicts the observation wells and the calculation done for the target period of 1963-2016 (53 years). When a fair fit between the observed and measured heads was satisfactorily reached, model calibration was terminated. Because of the small number and length of the observation results, the model calibration in this analysis can be considered based on these available data. When more evidence on data becomes usable in the future, further model calibration can be applied and the model results can be improved.  Hence, the model was calibrated to work under transient state conditions. The simulation in this case computes piezometric surface variations over time; hence, the parameters involved in the temporal equation can be calibrated (either a specific yield or a specific storage). Aquifer boundary inflows (flux at the boundary, meteoric recharge) and outflows (withdrawals) were also assessed and analyzed. The goal was to replicate the piezometric oscillations found at the control points as accurately as possible during the measurement period. The reference period for transient simulations corresponds to the one when GW level measurements were made. The computed piezometric heads found in steady-state regime acted as the initial conditions. The model was set up with the initial concentrations of SWI in 1961/1962. The simulation lasted a total of 20,075 days, beginning in 1962 and ending in 2016.
These simulations make use of wells that pump water mostly used for irrigation and DWS. Irrigation started in 1981 and lasted until 2016; withdrawals differ over time and are allocated based on both climate conditions and piezometric variations. The Thornthwaite technique was used to assess the monthly infiltration value that was distributed evenly throughout the month based on the meteoric recharge from the effective monthly precipitation. After calibration, the adopted basic storage was estimated to be between 0.1 and 4%. The calibration findings in a transient state regime are presented in the plots (Figure 17), which compare observed and computed heads. For the simulation time under consideration, this distinction demonstrates adequate consensus between measured and computed heads in various observation wells. GW pumping produces a significant drop in the GW level (piezometers 342/3, 438/3, and 1380/3), while irrigation in 1981 results in an increase in the water level (piezometers 32/3, 1432/3, and 1685/3).

Water Budget from 1961 to 2016
We assessed the amounts of the water budget at both the model borders and inside the aquifer itself. The simulation shows the start of the SWI and its progression over time, the pollution concentration, the intruded SWI length, and other mass balance elements.
The aquifer system has a negative balance of over 12 hm 3 (approximately 33,000 m 3 /day) between 1961/62 and 2015/16. It is worth noting that this time was marked by lowerthan-average meteoric recharge (−60%). The quantity of recharge (varies between 54 hm 3 and 20 hm 3 ) and the quantity of water crossing the alluvial plain (approximately 20 hm 3 inflow and over 25 hm 3 outflows) are the most interesting key components. The disparity (approximately 5 hm 3 ) reflects the net outflow of GW into the neighboring alluvial plain, which happens mostly on the left bank of the Loukkos wadi, as well as the amount of water outflowing towards the Atlantic Ocean (about 11 hm 3 ). Furthermore, the amount of GW pumping is approximately 20 hm 3 . Figure 18a depicts time differences in groudwater flow for each part of the water budget; note that substantial meteoric infiltration occurs primarily in 1962, 1968, 1976, 1996, and 2009 and that the piezometric level was highly affected by this recharge during these times.
The GW balance between 1962 and 2016 indicates that the aquifer had more freshwater storage between 1962 and 1976 and a decrease in aquifer recharge associated with EOP occurring between 1980 and 1990. As a result, seawater moved inland in 1982 and 1991 (Figure 18b,c). It also demonstrates that a decrease in recharge and EOP, especially in 1998, 2004, 2011, and 2015, increased SWI, though less pronounced than the intrusion of seawater in 1982 and 1991. The areas that have been intruded are as follows: (1) Between 1976 and 2016, the SWI entered the first region west of the ONEE pumping wells. Following that, the invaded region stretched along the coastal line for about 6 km in length and 0.5 km in width (see Figure 18, the small orange frame). (2) The aquifer is contaminated in the NW coastal plain, where the toe extends some 0.5 km inland. The contamination of the aquifer is limited beyond these areas.
Hence, when seawater directly joins the aquifer environment, it allows the chemical composition of GW to deteriorate by SWI. The situation is exacerbated by the presence of seawater at some times of the year when the piezometric levels are lower. The SWI edge extended 0.5 km into the aquifer bottom. We also observe that SWI increased in the northwestern sector in 1991/92, owing primarily to (1) decreased recharge induced by CC and intermittent droughts and (2) overexploitation of GW by extensive pumping out from the aquifer system for the water supply of the city of Larache, rural areas, and irrigation.

Climate Change, Over-Pumping, and SLR Impacts
The density-dependent numerical model that was designed was completed to simulate the flow and transport from 2017 to 2040, based on climate projections and GW management scenarios established by the National Office of Drinking Water (ONEE) in 2016 and ranging from 2017 to 2040 to meet the future water needs of both urban centers, Larache and Ksar El Kebir cities. Climate projections, under RCP 4.5, indicate a temperature rise of about 0.45 • C and a 16.7% decrease in precipitation from 2016 to 2050. Furthermore, the sea level will grow from 7 cm in 2020 to 15 cm by 2040.
For a period of 24 years, three planning scenario schemes were employed to simulate future changes in drawdown and salinity concentrations ( Figure 19). The first scenario assumes that the same conditions are maintained and that the aquifer pumping rate of approximately 21.52 hm 3 /year is maintained until 2040. Surface water or a desalination plant will meet the future increased water demand. The change in the aquifer's GW quality is also examined in order to determine the area affected by SWI. The second scenario involves increasing pumping rates until 2040 to supply the growing water demand of the Larache population. The third scenario assumes that additional GW abstractions will be required until 2040 to supply the water demands of both urban centers, Larache and Ksar El Kebir. Figure 19 depicts the three pumping scenarios and their progression to 2040. The projected drawdown and seawater volumes intruding into the aquifer are depicted in Figure 20, Figure 21 and Table 6.   Scenario 1, which assumes that the current pumping is maintained and future water demand will be provided by surface water or by a desalination plant, has less influence on the renewable resources and the water quality of the ROOCA. Indeed, from 2020 to 2040, we will note a drop in SWI volume (Figure 20b), which is directly related to an increase in hydraulic head (Figure 20a) due to increased predicted recharge from 2017 to 2040 in the study area. Figure 20c also demonstrates that the salinity concentration will be almost zero at piezometer 1380/3. Scenario 2 depicts a significant increase in salinity in the northwestern area, closer to the shoreline. The maximum extent of SWI will increase to 3.5 km deeper in the aquifer. The SWI volume intruding into the aquifer continues to rise (Figure 20b), while the GW level continues to decline, reaching maximum values around -10 m (Figure 20a). In 2040, however, the seawater would not reach the first series of wells. At piezometer 1380/3, the salinity concentration is also projected to rise to reach values around 20 g/L (Figure 20c). Scenario 3 is the pessimist one and also shows the predicted drawdown and saltwater extent in the aquifer by 2040 depicted in Figure 21. Note that sea salt concentration (35 g/L) is in red and freshwater concentration, almost 0 g/L, is in blue. There will be a greater decrease in hydraulic heads because of intensive pumping discharge to address the two cities' water needs. Indeed, hydraulic heads will hit negative values of around −20 m in the main ONEE well field sector (Figure 20a), where the drawdown will rise by 25 m, which is a substantial increase. The aquifer is also contaminated by SWI in the coastal part's northwest region, where the toe would reach about 5.2 km inland with an invaded area of about 31 km 2 ( Figure 21a) and would reach high salinity (15-25 g/L) in Layer 3 (Figure 21b). As a consequence, 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 2 g/L, will be already reached in 2032 and will continue to increase up to 6 g/L in 2040 (Figure 20a). A comparison between the second and third scenarios indicates that the salinity is expected to rise by 2040 from 20 to 30 g/L at piezometer 1380/3 located 1.5 km away from the coast (Figure 20c).

Conclusions
GW is the important source of freshwater in the Rmel-Oulad Ogbane coastal plain in the low Loukkos basin in Morocco, where 98% of the domestic water, as well as the entire industrial water supply, is dependent on GW. Therefore, it is imperative to protect GW from SWI in this area. Hence, there is a need for tools that can guide and assist the manager in decision-making regarding the use, management, and planning of water resources. Currently, these decision support elements are provided by efficient technical tools such as GIS, geostatistical analysis, and conceptual and mathematical models, as developed in this research.
Using data from boreholes and hydrogeological investigations, we modeled the aquifer reservoir using a 3D GSIS. Then, a geostatistical model was produced from physicochemical, piezometric, and hydrodynamic parameters. These outputs were used to update the water balance in 2013/2014 and to develop a good conceptual model of the aquifer. Finally, to predict the current extent of SWI and provide useful information for the protection of GW resources, a three-dimensional numerical model of density-dependent GW flow and miscible salt transport of the subsurface aquifer was developed to assess the current extent of SWI in the study area with the aim of determining the impacts of CC due to increasing of temperatures, decreasing precipitation, and SLR during the 21st century. The developed model incorporated regional geologic, geographic, and hydrogeological features. The model input parameters were determined from analysis of well logs, well driller's reports, and pumping tests. All these inputs were used to simulate three-dimensional variable-density GW flow under steady and transient states.
Due to the scarcity of some observed water quality data and continuously monitored head data, more field measurements, such as vertical salinity profile and trace element studies for dispersivity, need to be performed in order to improve the reliability of the model. For this performed calibration, a total of 22 observed hydraulic head values were used, resulting in good agreements between the observed and calculated hydraulic heads. When more data become available in the future, additional calibration would be needed. Moreover, an optimization model for rational management of the aquifer must be developed.
Climate projections used in this assessment comprise a combination of regional climate modeling projection data, generated from RICCAR, and a set of local observation datasets for precipitation and temperature for the study area. Projected temperature and precipitation were obtained from CNRM-CM5, EC-EARTH, and GFDL-ESM2M RCMs under RCP 4.5 and RCP 8.5 scenarios. These projections show a decrease in precipitation and an increase in temperature for both scenarios. As a result, CC would almost certainly have a detrimental effect on SLR, reducing the availability of new GW resources. The increase in sea level would affect coastal aquifers, shifting the saltwater interface farther inland. Indeed, the model took into account the predicted SLR and used it to adjust the boundary conditions. The variation in recharge was determined by taking into account the variations in return from irrigation and climatic parameters (precipitation and temperature). The numerical simulations were conducted for a period of approximately 76 years and dealt with SWI relating to GW abstraction, climatic parameters, and SLR. The simulation results under RCP 4.5 show that the maximum extent (about 5.2 km) of SWI would increase in 2040 in the northwestern sector of the study area. The water quality would be most affected in the ONEE pumping area, which is directly adjacent to the seashore. As a result, the GW abstraction associated with CC is the primary driver of SWI in the study area. Furthermore, the reduction in recharge and the rise in sea level caused by CC exacerbate saltwater intrusion into the aquifers, reducing the fresh GW resources.
The primary impact of this SWI in the ROOCA would be unnecessary over-pumping that would deplete renewable water resources. However, this situation can be improved by the use of surface water for irrigation (provided from the neighboring dam reservoir), a desalination plant project for DWS, and artificial recharge of the aquifer. GW recharge with recycled water would be also an effective and feasible way to address the rapid GW depletion and saltwater intrusion in the ROOCA. Recycled water is a sustainable and reliable source of local water that should be viewed as a valuable resource. Indeed, GW recharge is an excellent utilization of recycled water as it provides natural storage (which allows for drought mitigation or withdrawal when demand for water increases) and soil treatment (with surface spreading), and it can be used to directly prevent SWI (with a direct injection barrier, for instance). Thus, the city of Larache, which is located directly at the northern limit of the study site, has a good potential for wastewater to be produced for the artificial recharge of the coastal fringe of the aquifer, which is distinguished by sand dunes, very favorable to infiltration and natural purification.
This would greatly increase the GW production in the coastal sectors of the aquifer and would protect freshwater from SWI. Such long-term results and findings will help the local decision-makers and all relevant stakeholders to better plan, manage, and improve the fresh GW resources for the ROOCA.
Morocco has more than 3500 km of coastline (two maritime facades: Atlantic, 2,934 km; Mediterranean, 512 km) with several thousand hectares of coastal plains, where irrigated agriculture is well developed, such as in this case study, in addition to more industries and the most important cities in the country (Casablanca, Rabat, Tangier, Agadir, etc.). These coastal plains contain coastal aquifers that are overexploited and threatened by SWI. Some large cities are already supplied with surface water, GW, and desalination water (Agadir, Casablanca, Al Hoceima). Therefore, this study will serve as a pilot study in order to implement the established methodology. It is also recommended to complement this study with an in-depth study on the choice of artificial recharge sites and a study to optimize the management of conventional and unconventional water resources in order to minimize the energy costs associated with desalination and wastewater treatment, especially during wet and dry periods caused by CC.
It is also recommended to extend this established methodology to study further similar coastal aquifers abroad in terms of climate conditions, such as in the Mediterranean region, as this will help the decision-makers in water resources planning and development and securing sustainable GW management of the coastal aquifers. This methodology can be completed by an impact modeling study based on the application of artificial recharge (from available sources) in well-chosen sites in order to improve the production of water resources and limit the entry of seawater into freshwater of the coastal aquifer. Acknowledgments: This study is a part of the first author's Ph.D. thesis at the Mohammadia School of Engineers, UM5R, as well as a component of a research collaboration between the Moroccan Ministry of Water and the regional water center of Maghreb at EMI. We also thank DRPE and ABHL (Loukkos) for the remarkable cooperation that has occurred between both the university and the administrations and for providing data that have been used in this study. Partial financial support was also provided by the FOCP of Casablanca.

Conflicts of Interest:
The authors declare no conflict of interest.